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1 Introduction 

1.1 1925: Einstein's prediction for the ideal Bose gas 

Einstein considered N non-interacting bosonic and non-relativistic particles in a cubic 
box of volume with periodic boundary conditions. In the thermodynamic limit, defined 
as 

N 

iV, L ^ oo with — = p = constant, (1) 
a phase transition occurs at a temperature Tc defined by: 

pAL(Tc) = C(3/2) = 2.612... (2) 

where we have defined the thermal de Broglie wavelength of the gas as function of the 
temperature T : 

and where C{a) = Z^fcLi 1/^" is the Riemann Zeta function. 

The order parameter of this phase transition is the fraction Nq/N of particles in 
the ground state of the box, that is in the plane wave with momentum p = . For 
temperatures lower than this fraction Nq/N remains finite at the thermodynamic 
limit, whereas it tends to zero when T > Tc : 

T>Tc ^-0 (4) 

^ — (^) 

For T < Tc the system has formed a Bose- Einstein condensate in p = . The number 
Nq of particles in the condensate is on the order of N , that is macroscopic. As we will 
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see, the macroscopic population of a single quantum state is the key feature of a Bose- 
Einstein condensate, and gives rise to interesting properties, e.g. coherence (as for the 
laser) . 

1.2 Experimental proof? 

The major problem encountered experimentally to verify Einstein's predictions is that at 
densities and temperatures required by Eq.(|D at thermodynamic equilibrium almost all 
materials are in the solid state. 

An exception is He ^ which is a fluid at T = . However He ^ is a strongly interacting 
system. In He ^ in sharp contrast with the prediction for the ideal gas Eq.@, Nq/N < 
10% even at zero temperature Q 

The solution which victoriously led to Bose-Einstein condensation in atomic gases is 
to bring the system to extremely low densities (much lower than in a normal gas) and to 
cool it rapidly enough so that it has no time to recombine and solidify. The price to pay 
for an ultralow density is the necessity to cool at extremely low temperatures. Typically 
one has in the experiments with condensates: 

p < lO^^atoms/cm^ (6) 
T < l/jK. (7) 

The critical temperatures range from 20 nK to the /i K range. 

Bose-Einstein condensation was achieved for the first time in atomic gases in 1995. 

The group of Eric Cornell and Carl Wieman at JILA was first, with Rb atoms 0. 

They were closely followed by the group of Wolfgang Ketterle at MIT with Na atoms 

and the group of Randy Hulet at Rice University with Li atoms . Nowadays there 

are many condensates mainly with rubidium or sodium atoms. No other alkali atoms 

than the ones of year 1995 has been condensed. Atomic hydrogen has been condensed in 

1998 at MIT in the group of Dan Kleppner 0; the experiments on hydrogen were actually 

the first ones to start and played a fundamental pioneering role in developing many of 

"^Amusingly the ideal gas prediction Eq.(H) does not give a too wrong result for the transition temper- 
ature in helium. Note that the condensate fraction Nq/N should not be confused with the superfluid 
fraction: at T = the superfluid fraction is equal to unity. 
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the experimental techniques having led the alkalis atoms to success, such as magnetic 
trapping and evaporative cooling of atoms. 

In our lectures we do not consider the experimental techniques used to obtained and to 
study Bose-Einstein condensates as they are treated in the lectures of Wolgang Ketterle 
and of Eric Cornell at this school. 

1.3 Why interesting? 

1.3.1 Simple systems for the theory 

An important theoretical frame for Bose-Einstein condensation in interacting systems 
was developed in the 50's by Behaev, Bogohubov, Gross, Pitaevskii in the context of 
superfluid helium. This theory however is supposed to work better if applied to Bose 
condensed gases where the interactions are much weaker. 

The interactions in ultracold atomic gases can be described by a single parameter 
a , the so-called scattering length, as interactions take place between atoms with very 
low relative kinetic energy. The gaseous condensates are dilute systems as the mean 
interparticle separation is much larger than the scattering length a : 

p\a\^ < 1. (8) 

This provides a small parameter to the theory and, as we shall see, simple mean field 
approaches can be used with success to describe most of the properties of the atomic 
condensates. 

1.3.2 New features 

Atomic gases offer some new interesting features with respect to superffuid helium 4: 

• Spatial inhomogeneity. This feature can be used as a tool to detect the presence of 
a Bose-Einstein condensate inside the trap: in an inhomogeneous gas Bose-Einstein 
condensation occurs not only in momentum space but also in position space! 

• Finite size effects: The number of atoms in condensates of alkali gases is usually 
A^o < 10^ . The hydrogen condensate obtained at MIT by Kleppner is larger 
Nq ~ 10^ . Interesting finite size effects, that is effects which disappear at the 
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thermodynamic limit, such as Bose-Einstein condensates with effective attractive 
interactions ( a < ), can be studied in relatively small condensates. 

It is also interesting to consider small condensates where some interesting quan- 
tum aspects concerning coherence properties of the condensates, such as collapses 
and revivals of the relative phase between two condensates [^, could perhaps be 
measured [0]. 

• Tunability: Condensates in atomic gases can be manipulated and studied using the 
powerful techniques of atomic physics (see the lectures of Wolfgang Ketterle and 
Eric Cornell). Almost all the parameters can be controlled at will, including the 
interaction strength a between the particles. The atoms can be imaged not only 
in position space, but also in momentum space, allowing one to see the momentum 
distribution of atoms in the condensate! One can also tailor the shape and intensity 
of the trapping potential containing the condensate. 



2 The ideal Bose gas in a trap 

Let us consider a gas of non-interacting bosonic particles trapped in a potential f/(r) 
at thermal equilibrium. As the particles do not interact thermal equilibrium has to be 
provided by coupling to an external reservoir. In the grand-canonical ensemble the state 
of the gas is described by the equilibrium N -body density matrix 

p=^exp[-(3(H~fiN)] (9) 

where S is a normalization factor, H is the Hamiltonian containing the kinetic energy 
and trapping potential energy of all the particles, is the operator giving the total 
number of particles, (3 = 1/kBT where T is the temperature, and n is the chemical 
potential. One more conveniently introduces the fugacity: 

z = exp[(3fi]. (10) 

2.1 Bose-Einstein condensation in a harmonic trap 

Let us consider the case of a harmonic trapping potential U{f) : 

Uif) = \m{uy + uy + ulz^). (11) 
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We wish to determine the properties of the trapped gas at thermal equihbrium; the 
calculations can be done in the basis of harmonic levels or in position space. 

2.1.1 In the basis of harmonic levels 

Let us consider the single particle eigenstates of the harmonic potential with eigenvalues 
labeled by the vector: 

— * 

l^{lx,ly,lz) la^0,l,2,3... {a^ x,y,z). (12) 

One has: 

ef— IxhoJx + lyflOJy + Iz^Z (13) 

where the zero-point energy {h/2){ujx -\- ujy -\- uj^) has been absorbed for convenience in 
the definition of the chemical potential. Let us consider the case of an isotropic potential 
for which all the uja 's are equal to a; , so that e^— Ifiuj with I = + ly + Iz ■ 

The mean occupation number of each single particle eigenstate in the trap is given by 
the Bose distribution: 

■1 



' exp[/3(e^--//)] 



z 



■ exp(/3lhuj) — 1 



(14) 



Since n^- has to remain positive (for Z =0,1,2 ...), the range of variation of the fugacity 
z is given by 

0<z<l. (15) 

The average total number of particles N is obtained by summing over all the occupation 
numbers: N = I]f^f > ^ relation that can be used in principle to eliminate z in terms 
of . It is useful to keep in mind that for a fixed temperature T , N is an increasing 
function of z . 

In the limit z — one recovers Boltzmann statistics: n^- oc exp(— /Se^-) . We are 
interested here in the opposite, quantum degenerate limit where the occupation number 
of the ground state Z = of the trap, given by 

z 



No^n^^-—, (16) 

J- z 



diverges when z ^ 1 , which indicates the presence of a Bose- Einstein condensate in the 
ground state of the trap. 
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We wish to watch the formation of the condensate when z is getting closer to one, 
that is when one gradually increases the total number of particles N . The essence of 
Bose-Einstein condensation is actually the phenomenon of saturation of the population 
of the excited levels in the trap, a direct consequence of the Bose distribution function. 
Consider indeed the sum of the occupation numbers of the single particle excited states 
in the trap: 

N'-T^nr- (17) 
The key point is that for a given temperature T , N' is bounded from above: 

^' = E f- exp(/3mu;) - ll ' <Y1 [exp(/3ma;) - 1]"^ = . (18) 



z 



1^0 



Note that we can safely set z = 1 since the above sum excludes the term I — . 

If the temperature T is fixed and we start adding particles to the system, particles 
will be forced to pile up in the ground state of the trap when N > , where they will 
form a condensate. Let us now estimate the "critical" value of particle number A^^ax ■ 

We will restrict to the interesting regime ksT ^ Tiw : in this regime Bose statistics 
allows one to accumulate most of the particles in a single quantum state of the trap 
while having the system in contact with a thermostat at a temperature much higher 
than the quantum of oscillation hu . a very counter-intuitive result for someone used 
to Boltzmann statistics! On the contrary the regime ksT -C TiuJ would lead to a large 
occupation number of the ground state of the trap even for Boltzmann statistics. 

A first way to calculate N[^^^ is to realize that the generic term of the sum varies 
slowly with / as ksT ^ hoo so that one can replace the discrete sum I]f_^g by an 
integral /^^>o d^l . As we are in the case of a three-dimensional harmonic trap there is no 
divergence of the integral around / = . 

We will rather use a second method, which allows one to calculate also the first cor- 
rection to the leading term in UbT /fiuj . We use the series expansion 

-< — T OO 

-x-Y.^-''^ (19) 



1 1 — 



which leads to the following expression for A^4iax ) if one exchanges the summations over 
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/ and k : 



NL. = E E eM-muj E = E 

k=l a k=l 



1 — cxp[—Pfiujk] 



- 1 



(20) 



We now expand the expression inside the brackets for small x 

3 

I -1 



1 — exp[— x] 
and we sum term by term to obtain 



1 3 



2x' 



+ ... 



(21) 



NL 



huj , 



C(3) + ^ 



3 (kuT" 



C(2) 



(22) 



Note that the exchange of summation over k and summation over the order of expansion 
in Eq. (^I]) is no longer allowed for the next term 1/x , which would lead to a logarithmic 
divergence (that one can cut "by hand" at A; ~ ksT/huj ). 

One then finds to leading order for the fraction of population in the single particle 
ground state: 



No N 



N' 

max 



1-C(3) 



knT' 



N N ' \huj 

where the critical temperature T° is defined by: 



1 

N 



1 



C(3) 



N 



(23) 



(24) 



and C(3) = 1.202... . Note that the universal law ([23| ) differs from the one obtained in 
the homogeneous case (H) usually considered in the literature. 

The present calculation is easily extended to the case of an anisotropic harmonic trap. 
To leading order one finds 



max 



ksT 

flUJ 



C(3) 



(25) 



where uj 



n/3 



is the geometric mean of the trap frequencies. One can also 



calculate A^4iax two-dimensional and one-dimensional models. One also finds in these 
cases a finite value for A^4iax • ^^e saturation of population in the single particle excited 
states applies as well and one can form a condensate, a situation very different from the 
thermodynamical limit in the homogeneous ID and 2D cases. 
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2.1.2 Comparison with the exact calculation 



One can see in figure |T] that the first two terms in the expansion Eg . (^21) , combined with 
the approximation Ai"o/A^ — 1 — N'^.-^^/N , give a very good approximation to the exact 
condensate fraction for = 1000 particles only. 



o 




0.4 0.6 O.B 

T/T,° 



Figure 1 : Condensate fraction versus temperature for an ideal Bose gas in a spherically symmet- 
ric trap with A'^ = 1000 particles. The circles correspond to the exact quantum calculation. The 
solid line corresponds to the prediction Nq/N 0:^ with A^^^x given by the two terms 

in the expansion Eq.(^). The dashed line corresponds to the prediction Nq/N ~ 1 — N^^^/N 
given by the leading term in Eq.(^). This figure was taken from 



with A^^^^ 



2.1.3 In position space 

A very important object in the description of the state of the gas is the so-called one-body 
density matrix. We can define it as follows. 
Consider a one-body observable 



N 

1=1 



(26) 
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where X{i) is the observable for particle number i and where is the operator giv- 
ing the total number of particles. The one-body density matrix pi is defined by the 
requirement that for any X : 

{X) ^ Tr[piX(l)]. (27) 

For the particular case of X equal to the identity it follows X = N and {X) = Tr[pi] = 
(A^) so that our one-body density matrix is normalized to the mean number of particles 
in the system. 

An equivalent definition of pi in the second quantized formalism is simply 

{f'\pi\f) = {ij\r)ip{r')) (28) 

where ilj{f) is the atomic field operator, annihilating an atom in r . 

At thermal equilibrium in the grand-canonical ensemble, the one-body density matrix 
of the ideal Bose gas is given by 

Pi = (29) 

z ^ exp{i3hi) — 1 

where the single-particle Hamiltonian in the case of a spherically symmetric harmonic 
trap is 

hi = — 1 — mcj^r ^ fiuj. (30) 

2m 2 2 

Here again we have subtracted the zero-point energy for convenience. The Bose formula 
Eq.(p!^ corresponds to the diagonal element of pi in the eigenbasis of the harmonic 
oscillator (the off-diagonal elements of course vanish). In position space the diagonal 
term 

(f |pi|f) = p(f) (31) 

gives the mean spatial density of the gas. 

In order to calculate the density we use the series expansion Eq.(p!9D to rewrite pi as 
follows: 



oo 



k=l 

This writing takes advantage of the fact that the matrix elements 

(f|e-^'='^i|f') (33) 
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are known for an harmonic oscillator potential 



p(f; 



k=l 



[1 — exp{—2f3khuj)) '^^^ exp 



One then obtains explicitly 

(3khuj\ 



muor 



h 



tanh 



(34) 



One can identify the contribution of the condensate to this sum when z — 1^ . When 
the summation index k is large, what determines the convergence of the series is indeed 
the factor . Replacing the other factors in the summand by their asymptotic value for 
k +00 we identify the diverging part when z = 1 : 



A;=l 



moor 



1 -z 



|0O,O,o(Ol' = ^o|0O,O,o(OP 



(35) 



where 0o,o,o('^) is the ground state wave function of the harmonic oscillator. 

Numerically we have calculated the total density p(r ) for a fixed temperature ksT = 
20huj and for increasing number of particles (see figure 0). Here the maximal number of 
particles one can put in the excited states of the trap is A^^ax — C(3)(^s^/^'-^)^ — 10^ • 
When -C N^^^ the effect of an increase of is mainly to multiply the density 
by some global factor (the curves in logarithmic scale in figure |^ are parallel one to 
the other). When is becoming larger than A^4iax ^ peak in density grows around 
r = , indicating the formation of the condensate, whereas the far wings of the density 
distribution saturate, which reflects the saturation of the population of the excited levels 
of the trap. 



2.1.4 Relation to Einstein's condition pX^B = C(3/2) 

In the limit ksT ^ hu we can actually calculate the value Pmax(^) which the density 
p'(r ) of particles in the excited states of the trap saturates when z — 1 . We simply use 
the expansion Eq. (|3^) , subtracting from the total density p(r ) the contribution of the 
condensate A'o|0o,o,o(^)P • The resulting series is converging even for z = 1 so that we 
can take safely the semiclassical limit ksT ^ hu term by term in the sum: 

1 oo ^fc 



p'{f) 



E 



ti k'/' 



exp 



--kPmu'^r^ 



X3 



93/2 



;exp ^— /J-mcu^r^ 



(36) 



where 
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10 20 



r[ao] 

Figure 2: Spatial density for an ideal Bose gas at thermal equilibrium in a harmonic trap of 
frequency uj . The temperature is fixed to ksT = 20huj and the number of particles ranges 
from = 500 to = 32000 between the lowest curve and the upper curve, with a geometrical 
reason equal to 2. The unit of length for the figure is oq = {fi/2mujY^'^ , that is the spatial 
radius of the ground state of the trap. 

We term this approximation semiclassical as (i) one can imagine that the classical limit 
— > is taken in each term k of the sum, giving the usual Gaussian distribution for the 
density of a classical harmonic oscillator at temperature kBT/k , but (ii) the distribution 
still reflects the quantum Bose statistics. 

If now we set z = 1 in (|36D to express the fact that a condensate is formed we obtain 

PLAr = 0) ^ -^93/2(1) = (38) 

We therefore recover Einstein's condition provided one replaces the density p of the 
homogeneous case by the density at the center of the trap. 

2.2 Bose-Einstein condensation in a more general trap 

We now extend the idea of the previous semiclassical limit to more general non-harmonic 
potentials. This allows to find the condition for Bose-Einstein condensation in presence 
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of a non-harmonic potential. This will prove useful in presence of interactions between 
the particles where the non-harmonicity is provided by the mean field potential. 



2.2.1 The Wigner distribution 

The idea is to find a representation of the one-body density matrix having a simple (non 
pathological) behavior when h ^ . Let us take as an example a single harmonic 
oscillator. The density matrix is then of the form: 

1 



a = —e 'f^^ho 
Z 



(39) 



where H^^o is the harmonic oscillator Hamiltonian. As shown in all the matrix elements 
of a can be calculated exactly: 



{f\a\f') 



(27r)3/2(Ar)3 



exp 



2(Ar)2 



exp 



(r — r ' )^ 



2e 



The relevant length scales are the spatial width of the cloud Ar 

(Ar)2 

and the coherence length ^ : 



-cotanh 



2muj 



2kBT 



2 2h I huj 

4 = tanh 



muj 



2kBT 



(40) 



(41) 



(42) 



If we now take the classical limit ^ ^ (in more physical terms the limit hu ^ ksT 
then: 



(Ar)2 



X2 



(43) 
(44) 



rnksT 27r 

In the limit — > the h dependence of causes {f\a\f')^0 for fixed values of r, r' 
unless r = r' : the limit is singular. 

To avoid this problem one can use the Wigner representation of the density matrix, 
introduced also in the lectures of Zurek and Paz: 

(Pu u 



(r 



cr 



2 



(45) 
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The Wigner distribution is the quantum analog of the classical phase space distribution. 
In particular one can check that the Wigner distribution is normalized to unity and that 



d^rW{r,p) 
(f'pW{r,p) 



{p\a\p) 
{f\a\f). 



(46) 
(47) 



An important caveat is that W is not necessarily positive. 

For the harmonic oscillator at thermal equilibrium the integral over u in Eq. (^5]) is 
Gaussian and can be performed exactly: 



W{f,p) 



exp( 



(27rArAp)3 2(Ar)2 
where Ap = h/^ . If we take now the limit /i : 



)exp( 



p 



2{Apy 



(Ar)2 
(Apr 



raksT 

so that Wif^p) tends to the classical phase space density. 



(48) 

(49) 
(50) 



2.2.2 Critical temperature in the semiclassical limit 

Let us turn back to our problem of trapped atoms in a non-harmonic trap where the 
single particle Hamiltonian is given by 

(51) 



p' 



and the one-body density matrix is given by Eq. (^2]) . For /i — > we have: 



W\e 



-^'^'^^](f,p-)^-exp 



(52) 



As we did before we put apart the contribution of the condensate. One then gets for 
the one-body density matrix of the non-condensed fraction of the gas in the semiclassical 
limit: 

1 +°° r / 



1 +CO 

IsL^ exp 

k=l 



p^ 



z 



exp 



p^ 



(53) 
(54) 
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We are now interested in the spatial density of the non-condensed particles in the 



semiclassical limit. By integrating Eg . (1531) over p we obtain: 



P'scin = y^^73/2(ze-^^('^)) (55) 

where Qa is defined in Eq. (|37D . The condition for Bose-Einstein condensation is z — *^ 
e^'^min where f/min = TamffU{f) is the minimal value of the trapping potential, achieved 
in the point f^^i^ . For z = e^^min the semiclassical approximation for the non-condensed 
density gives in this point: 

PL(^min) = T^^3/2(l) (56) 



^dB 



or 



pXls = 2.612... (57) 



Again Einstein's formula is recovered with p being the maximal density of the non- 
condensed cloud, that is the non-condensed density at the center of the trap. 

The semiclassical calculation that we have just presented was initially put forward 
in |]10|. We do not discuss in details the validity of this semi-classical approximation. 
Intuitively a necessary condition is ksT ^ AE where AE is the maximal level spacing 
of the single particle Hamiltonian among the states thermally populated. Some situations, 
where the trapping potential is not just a single well, may actually require more care. The 
case of Bose-Einstein condensation in a periodic potential is an interesting example that 
we leave as an exercise to the reader. 



2.3 Is the ideal Bose gas model sufficient: experimental verdict 
2.3.1 Condensed fraction as a function of temperature 

The groups at MIT and JILA have measured the condensate fraction Nq/N as function 
of temperature for a typical number of particles = 10^ or larger. We reproduce 
here the results of JILA (see figure This figure shows that the leading order 
prediction of the ideal Bose gas Eq.(^3|) is quite good, even if there is a clear indication 
from the experimental data that the actual transition temperature is lower than T° . This 
deviation may be due to finite size effects and interaction effects but the large experimental 
error has not allowed yet a fully quantitative comparison to theory. 
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Figure 3: Condensate fraction Nq/N as function of T/T^ where is the leading order 
ideal Bose gas prediction Eq.(p^. Circles are the experimental results of |11] while the dashed 
line is Eq.(|3|). 



2.3.2 Energy of the gas as function of temperature and number of particles 

In the experiments one produces first a Bose condensed gas at thermal equilibrium. Then 
one switches off suddenly the trapping potential. The cloud then expands ballistically, 
and after a time long enough that the expansion velocity has reached a steady state value 
one measures the kinetic energy of the expanding cloud. 

Suppose that the trap is switched off at t = . For t = 0~ the total energy of the 
gas can be written as 

-£^101(0") = -Ekin + -Satrap + -E'int, (58) 

that is as the sum of kinetic energy, trapping potential energy and interaction energy. At 
time t = 0'^ there is no trapping potential anymore so that the total energy of the gas 
reduces to 

Etot(0+) = ^kin + ^int. (59) 

In the limit t +00 the gas expands, the density and therefore the interaction energy 
drop, and all the energy Etot{0~^) is converted into kinetic energy, which is measured. 
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In figure § we show the results of JILA for Etot{0^) for temperatures around T° 
[1T| together with the ideal Bose gas prediction. The main feature of the ideal Bose gas 
prediction is a change in the slope of the energy as function of temperature when T 
crosses . One observes indeed a change of slope in the experimental results (see the 
magnified inset)! 

For T > Tf. the ideal Bose gas model is in good agreement with the experiment. For 
T < Tc we observe however that the experiment significantly deviates from the ideal Bose 
gas. 




Figure 4: Expansion energy of the gas Etot{0'^) per particle and in units of ksT^ as function 



of the temperature in units of . The disks correspond to the experimental results of |11]. 
The straight sohd hne is the prediction of Boltzmann statistics. The dashed curve exhibiting a 
change of slope is the ideal Bose gas prediction. The curved solid line is a piecewise polynomial 
fit to the data. The inset is a magnification showing the change of slope of the energy as function 
of T close to T = . The figure is taken from 



What happens at even lower values of T /T^ ? We show in figure ^ the expansion 



energy of the condensate per particle in the regime of an almost pure condensate |T2 
This energy then depends almost only on the number of condensate particles Nq , m. a. 
non-linear fashion. This is in complete violation with the ideal Bose gas model, which 
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predicts an energy per particle in the condensate independent of Nq . More precisely 
the ideal Bose gas prediction would be h{ujx + Uy + uJz)/4: where the Ua 's are the trap 
frequencies. In units of fc^ this would be in the 10 nK range, an order of magnitude 
smaller than the measured values. 



200 - 




Number of Condensed Atoms Nq (10^) 

Figure 5: Expansion energy of the condensate per particle in the condensate, divided by ks , 
as a function of the number of particles in the condensate. The experiment is performed at 
temperatures T <^ Tc . The triangles correspond to cases where the non-condensed cloud was 
not visible experimentally. The disks correspond to case where the non-condensed cloud could 
be seen. The figure is taken from [^]. The solid line is a fit of the interacting Bose gas prediction 
of a 



2.3.3 Density profile of the condensate 

The group of Lene Hau at Rowland Institute has measured the density profile of the 
condensate in a cigar-shaped trap, along the weakly confining axis z of the trap. As 
imaging with a light beam is used the actual density obtained in the experiment is the 
density integrated along the direction y of propagation of the laser beam, plotted in 
figure M for x = as function of z [1^ . The measured profile is very different from and 
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much broader than the Gaussian density profile of the ground state wavefunction of the 
harmonic oscillator. 
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Figure 6: Column density profile (see text) of a condensate along the weak axis z of a cigar- 
shaped trap. The experimental results of (dots) are very different from the ideal Bose gas 
prediction (dashed line). The solid line corresponds to the theoretical prediction of §^. 



2.3.4 Response frequencies of the condensate 

By modulating the harmonic frequencies of the trapping potential one can excite breathing 
modes of the condensate. For example the group at MIT modulated the trap frequency 
along the slow axis z of a cigar-shaped trap and observed at T <^ Tc subsequent 
breathing of the condensate at a frequency 1.569(4)ti;z . This frequency is not an integer 
multiple of uJz and can therefore not be obtained in the ideal Bose gas model. 

In conclusion the ideal Bose gas model may be acceptable as long as no significant 
condensate has been formed. If a condensate is formed interaction effects become impor- 
tant, and dominant at T <^ Tc . This serves as a motivation to the next sections of this 
lecture, which will deal with the interacting Bose gas problem. 
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3 A model for the atomic interactions 

The previous section ^ has shown that the ideal Bose gas model is insufficient to explain 
the experimental results when a condensate is formed. In this section we choose the 
model potential to be used in this lecture to take into account the atomic interactions. 
The reader interested in a more careful discussion of real interaction potentials is referred 



to 14 . 



3.1 Reminder of scattering theory 

We consider two particles of mass m interacting in free space via the potential V{fl—f^) 
depending on the positions rl, r2 only through the relative vector rl — r2 . The center of 
mass of the two particles is then decoupled from their relative motion, and the evolution 
of the relative motion is governed by the Hamiltonian: 

Hrei = f- + Vir) (60) 

where r = ?^ — is the vector of coordinates of the relative motion, p = [pi — p2)/2 is 
the relative momentum and /i = m/2 is the reduced mass. We assume in what follows 
that the potential V{f) is vanishing in the limit r ^ oo . 

3.1.1 General results of scattering theory 

The scattering states iplf) of the relative motion of the two particles are the eigenstates 
of i^rci with positive energy E . Writing E = ^^fc^/2/i and multiplying the eigenvalue 
equation by 2fi/h^ we obtain 

{A + k'^)tP{r) = '^V{f)ip{r). (61) 

h 

One has also to specify boundary conditions on ip to get the full description of a scattering 
state. This is achieved by means of an integral formulation of the eigenvalue equation. 

• Integral equation 

To obtain the integral formulation of the scattering problem we write the right hand side 
of the eigenvalue equation Eq. (|6TD as a continuous sum of Dirac distributions: 

(A + k^)^{r) = f d^f''^V{f')^{f')6{f- f'). (62) 

J h 
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We then find a solution of tliis equation witli a single Dirac distribution on the right hand 
side: 

(A,-+A;2)^G(r) = 5(r-f') (63) 
having the form of an outgoing spherical wave for r ^ oo : 

]^ gifc|r-'r"l 

Mr) = - — r^^. (64) 

Att \r — r'\ 

This is actually a Green's function of the operator A + k"^ . The scattering state of the 
full problem can then be written as 

O,, r pik\f-f'\ 

m = Mr) - / d'f'——Vif')M-") ■ (65) 

Airn J \r — r \ 

The first term ipQ is the incoming free wave of the collision, solving (A + k'^)ipo = ; we 
simply assume here that the incoming wave is a plane wave of wavevector k : 

ipoir) = exp[ik ■ r ]. (66) 

The remaining part of ip is then simply the scattered wave. 

• Born expansion 

When the interaction potential is weak one sometimes expands the scattering state V 
in powers of . In the integral formulation Eq. (^) of the eigenvalue equation this 
corresponds to successive iterations of the integral, the approximation for ip at order 
n + 1 in V being obtained by replacing by its approximation at order n in the 
right-hand side of the integral equation. E.g. to zeroth order in V , ip = ipo , and to first 
order in V we get the so-called Born approximation: 

2/i /■ ,3^, e^'^l^^-'"'. 



V'Born(r ) = Mr) - / d'f' ^-^^V{f')Mr')- (67) 

A-Kh J \r ~r'\ 

3.1.2 Low energy limit for scattering by a finite range potential 

Some results can be obtained in a simple way when the potential V has a finite range 
b , that is when it vanishes when r > b . 

• asymptotic behavior for large r 
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As the integration over the variable r ' is hmited to a range of radius b one can expand 
the distance from r to r ' in powers of r when r ^ b : 



r-f'-n + 0{-) (68) 



where n = f/r is the direction of scattering. The neglected term, scaling as , has a 
negligible contribution to the phase exp[ik\f— f'\] when r ^ kb"^ . One then enters the 
asymptotic regime for ip : 

^(f ) = Mr) + ^fd^) + O (69) 

where the factor , the so-called scattering amplitude, does not depend on the distance 
r : 

If the mean distance between the particles in the gas, on the order of p^^^^ , where 
p is the density, lies in the asymptotic regime for ip (that is p^^/^ ^ b, kb"^ ) the effect 
of binary interactions on the macroscopic properties of the gas will be sensitive to the 
scattering amplitude , and no longer to the details of the scattering potential. This 
is the key property that we shall use later in this low density regime to replace the 
exact interaction potential by a model potential having approximately the same scattering 
amplitude. 

• limit of low energy collisions 

Another simplification comes from the fact that collisions take place at low energy in 
the Bose condensed g 3iSGSi clS is on the order of ksT in the thermal gas, k 

becomes small at low temperature. 

If kb the phase factor exp[— -i/cn-r '] becomes close to one in the integral Eq. (|70D 
giving the scattering amphtude. The scattering amplitude then no longer depends on 
the scattering direction n , the asymptotic part of the scattered wave becomes spherically 
symmetric (even if the scattering potential is not!): one then says that scattering takes 
place in the s -wave only. 

Going to the mathematical limit — > we get for the scattering amplitude: 

kin) ^ -a. (71) 



Model for interactions 



26 



The quantity a is the so-called scattering amplitude; it will be the only parameter of our 
theory describing the interactions between the particles, and our model potential will be 
adjusted to have the same scattering length as the exact potential. When k is going to 
zero, the scattering state converges to the zero energy scattering state, behaving for large 
r as 



A numerical calculation of this zero energy scattering state is an efficient way of calculating 
a for a given potential V . Note that there is of course no connection between a and 
b , except for particular potentials like the hard sphere potential. 

3.1.3 Power law potentials 

In real life the interaction potential between atoms is not of finite range, as it contains 
the Van der Waals tail scaling as for large r It is fortunately possible to show 

for the class of power- law potentials, scaling as 1/r" , that several of our conclusions, 
obtained in the finite range case, hold provided that n > 3 . E.g. in the limit of small 
k 's only the s -wave scattering survives, and has a well defined limit for — , 
allowing one to define the scattering length. 

3.2 The model potential used in this lecture 
3.2.1 Why not keep the exact interaction potential ? 

For alkali atoms the exact interaction potential has a repulsive hard core, is very deep (as 
deep as 10^ Kelvins times ks for ^^'^ Cs), has a minimum at a distance ri2 on the order 
of 6 A(for cesium), and contains many bound states corresponding to molecular states of 
two alkali atoms (see figure |^) . 

There are several disadvantages to use the exact interaction potential in a theoretical 
treatment of Bose-Einstein condensation: 

1. V is difficult to calculate precisely, and a small error on V may result in a large 
error on the scattering length a . In practice a is measured experimentally, and 
this is the most relevant information on V in the low density, low temperature 
limit. 

^ or even as if r is larger than the optical wavelength. 




(72) 
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Figure 7: Typical shape of the interaction potential between two atoms, as function of the 
interatomic distance ri2 . The numbers are indicative and correspond to cesium. 

2. the presence of bound states of V with a binding energy much smaller than the 
temperature of the gas (there are 9 orders of magnitude between the potential depth 
10^ K and the gas temperature ~ Ifj, K) clearly indicates that the Bose condensed 
gases are in a metastable state; at the experimental temperatures and densities 
the complete thermal equilibrium of the system would be a solid. Direct thermal 
equilibrium theory, such as the thermal -body density matrix exp[—i3H] , cannot 
therefore be used with V . This is why even in the exact Quantum Monte Carlo 
calculations performed for alkali gases V is replaced by a hard sphere potential. 
Such a complication was absent for liquid helium, where the well-known exact V 
can be used [|^. 

3. V can not be treated in the Born approximation, because it is very strongly re- 
pulsive at short distances and has many bound states: even if the scattering length 
was zero, one would have to resum the whole Born series to obtain the correct result 
[We recall that for a potential as gentle as a square well of radius b , the Born ap- 
proximation applies when the zero-point energy for confinement within a domain of 
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radius b , /i^/2/i6^ , is much larger than the potential depth, which implies that no 
bound state is present in the well.] As a consequence naive mean field approxima- 
tions, which neglect the correlations between particles due to interactions, implicitly 
relying on the Born approximation, cannot be used with the exact V . 

The key idea is therefore to replace the exact interaction potential by a model potential 
(i) having the same scattering properties at low energy, that is the same scattering length, 
and (ii) which should be treatable in the Born approximation, so that naive mean field 
approaches apply. 

The model potential satisfying these requirements with the minimal number of pa- 
rameters (one!) is the zero-range pseudo-potential initially introduced by Enrico Fermi 
T7|, |T3 and having the following action on any two-body wavefunction: 



(n, r^2\V\iJi^2) = gS{fl - ri) 



_dri2 

The factor g is the so-called coupling constant 



d 

(n2^1,2(n,?^"2)) 



(73) 

J ri2=0 



9 = a 74 

m 

where a is the scattering length of the exact potential. The pseudo-potential involves a 
Dirac distribution and a regularizing operator. 

• Effect of regularization 

When the wavefunction ■i/'i 2 is regular close to f[ = 7^2 , one can check that the regular- 
izing operator has no effect, so that the pseudo-potential can be viewed as a mere contact 
potential g6{fi — -rj.) . 

When the wavefunction '?/'i 2 has a l/ri2 divergence: 

Wi,2{ri,r2) = h regular (75) 

ri2 

where A is the function of the center of mass coordinates only the regularizing operator 
removes the diverging part: 



dri2 V n2 
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In this way we have extended the Hilbert space of the state vectors of the particles 
with wave functions diverging as l/r^ ; note that these wavefunctions remain square 
integrable, as the element of volume scales as 3D. As we shall see this l/ri2 

divergence is a consequence of the zero-range of the pseudo-potential. 



3.2.2 Scattering states of the pseudo-potential 



Turning back to the relative motion of two particles we now derive the scattering states of 
the pseudo-potential from the integral equation Eq. (|65D . As the pseudo-potential involves 
a Dirac S{f') the integral over r ' can be performed explicitly: 



_d_ 

df 



- (rV(r")) 



As the factor 



C 



^ ')) 



(77) 



(78) 



J r'=0 



does not depend on r we find that has the standard asymptotic behavior of a scat- 
tering state in r but everywhere in space, not only for large r . This is due to the 
zero-range of the pseudo-potential. To calculate C , we multiply Eq.([77D by r , we take 
the derivative with respect to r and set r to zero. On the left hand side we recover the 
constant C by definition. We finally obtain: 



C = 1- aCik 



(79) 



so that C = 1/(1 -|- ika) and the scattering states of the pseudo-potential are exactly 
given by 

^ „ „ikr 

(80) 



a 



^ ' 1 + ika r 

The corresponding scattering amplitude, 

a 



fk 



1 + ika 



(81) 



does not depend on the direction of scattering, so that the pseudo-potential scatters only 
in the s -wave, whatever the modulus k is. The scattering length of the pseudo-potential, 
—fk=o = ci , coincides with the one of the exact potential. 
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Finally we note that the total cross-section for scattering of identical bosons by the 
pseudo-potential is given by a Lorentzian in k , 

and that the pseudo-potential obeys the optical theorem. 
3.2.3 Bound states of the pseudo-potential 

As a mathematical curiosity we now point out that not only the scattering states but also 
the bound states of the pseudo-potential can be calculated. A first way of obtaining the 
bound states is a direct solution of Schrodinger's equation. A more amusing way is to use 
the following closure relation: 

(0I^A:)(l^d = l -Abound (83) 

where {ipj:) is the scattering state given in Eq.(pO[) and Pbound is the projector on the 
bound states of the pseudo-potential. 

In calculating the matrix elements of this closure relation between perfectly localized 
state vectors |r ) and \f') and using spherical coordinates for the integration over k 
one ultimately faces the following type of integrals: 

dk (84) 
-oo 1 -|- ika 

We calculate / using the residues formula, by extending the integration variable k to 
the complex plane and closing the contour of integration by a circle of infinite radius, 
which has to be in the upper half of the complex plane as r + r' > . As the integrand 
in / has a pole in k = i/a , we find that I vanishes for a < , as the pole is then in 
the lower half of the complex plane. For a > the pole gives a non-zero contribution to 
the integral: 

a 

Finally we find that Pbound = for a < , corresponding to the absence of bound 
states, and Pbound = I "Abound) (Abound I for a > 0, Corresponding to the existence of a 
single bound state: 

1 e-"'/" . . 
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From Schrodinger's equation, we find for tlie energy of the bound state: 

-Abound = ^- (87) 

The existence of a bound state for a > and its absence for a < is a paradoxical 
situation. As we shall see in the mean field approximation, the case a > corresponds to 
effective repulsive interactions between the atoms, whereas the case a < corresponds 
to effective attractive interactions. In the purely ID case, the situation is more intuitive, 
the potential qidS^x) having a bound state only in the effective attractive case giD < . 
This paradox in 3D comes from the non-intuitive effect of the regularizing operator (an 
operation not required in ID), which makes the pseudo-potential different from a delta 
potential; actually one can shown in 3D that a delta potential viewed as a limit of square 
well potentials with decreasing width b and constant area does not scattered in the limit 
b^O. 



3.3 Perturbative vs non-perturbative regimes for the pseudo- 
potential 

3.3.1 Regime of the Born approximation 

As we will use mean field approximations requiring that the scattering potential is treat- 
able in the Born approximation, we identify the regime of validity of the Born approxi- 
mation for the pseudo-potential. 

As we have seen in the previous subsection the integral equation for the scattering 
states of the pseudo-potential can be reduced to the equation for C : 

C = 1- tkaC, (88) 

the scattering state being given by 

_^ pikr 

^^(f) = e^'^-^^-aC— . (89) 

The Born expansion will then reduces to iterations of Eg . (|88|) . To zeroth order in the 
interaction potential, we obtain Cq = so that ip^ reduces to the incoming wave. To 
first order, we get the Born approximation 

Ci = l- ikaCo = 1. (90) 
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To second order and third order we obtain 

C2 = 1- ikaCi = l-ika (91) 
C3 = 1 - ikaC2 = l-ika + {ikaf (92) 

so that the Born expansion is a geometrical series expansion of the exact result C = 
1/(1 + ika) in powers of ika . 

The validity condition of the Born approximation is that the first order result is a 
small correction to the zeroth order result. For the scattering amplitude this requires 

k\a\ < 1. (93) 

For the scattering state this requires 

r > a. (94) 

If one takes for r the typical distance p~^^^ between the particles in the gas, where p 
is the density, this leads to 

p^/^|a| < 1. (95) 

• Are the conditions for the Born approximation satisfied in the experiments ? 

To estimate the order of magnitude of k we average fc^ over a Maxwell-Boltzmann 
distribution of atoms with a temperature T = l/x K, typically larger than the critical 
temperature for alkali gases; the average gives a root mean square for k equal to 

For Na atoms used at MIT, with a scattering length of 50 aeohr , where the Bohr radius 
is ttBohr = 0.53 A, we obtain Aka = 2 x 10~^ . For rubidium Rb atoms used at JILA, 
with a scattering length of llOaBohr , "we obtain Aka = 0.1 . 

In the case of an almost pure condensate in a trap, the typical k is given by the 
inverse of the size R of the condensate, as the condensate wavefunction is not very far 
from a minimum uncertainty state. Generally this results in a much smaller Ak than 
Eq. (P^D , as R is much larger than the thermal de Broglie wavelength. One could however 
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imagine a condensate in a very strongly confining trap, such that R would become close 
to a ; in this case, not yet realized, the mean field theory has to be revisited. 

We turn to the second condition Eg. (pSj) . The typical densities of condensates are on 
the order of 2 x 10^^ atoms per cm ^ . For the scattering length of sodium this leads 
to p^/'^a ~ 0.015 <^ 1 . For the scattering length of rubidium this leads to 
0.034 ^ 1 . Both conditions for the Born approximation applied to the pseudo-potential 
are therefore satisfied. 

3.3.2 Relevance of the pseudo-potential beyond the Born approximation 

Let us try to determine necessary validity conditions for the substitution of the exact 
interaction potential by the pseudo-potential. 

First one should be in a regime dominated by s -wave scattering, as the pseudo- 
potential neglects scattering in the other wave. This condition is easily satisfied in the 
fi K temperature range for Rb, Na. 

Second the scattering amplitude of the exact potential in s -wave should be well 
approximated by the pseudo-potential. For isotropic potentials vanishing for large r as 
l/r" , with n > 5 , the s -wave scattering amplitude has the following low k expansion: 

/r° = (97) 

where r^. is the so-called effective range of the potential. To this order in k the result of 
the pseudo-potential corresponds to the approximation re = . When r^. is on the order 
of a (which is the case for a hard sphere potential, but not necessarily true for a more 
general potential) the term in can be neglected if fc^rg <^ 1/a , that is (k.af' <^ 1 ; 
there is therefore no meaning to use the pseudo-potential beyond the Born regime. 

Consider now the case re -C |a| . The term r^k^ remains small as compared to 1/a 
for k\a\ < 1 . For k\a\ ^ 1 the term ik dominates over 1/a ; k'^r^ remains small as 
compared to ik as long as kr^ <^ 1 . The use of the pseudo-potential may then extend 
beyond the Born approximation. 

An example of a situation with re <^ |a| is the so-called zero energy resonance, where 
a is diverging. When a bound state of the interaction potential is arbitrarily close to the 
dissociation limit, the scattering length diverges a — > +oo , the bound state has a large 
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tail in r scaling as e ''/'^/r and the bound state energy scales as —h'^/ma'^ |T^, |20[ . 
These scaling laws hold for the pseudo-potential, as we have seen. 



4 Interacting Bose gas in the Hartree-Fock approxi- 
mation 

Now that we have identified a simple model interaction potential treatable in the Born 
approximation we use it in the simplest possible mean field approximation, the so-called 
Hartree-Fock approximation. This approximation was applied to trapped gases for the 
first time in 1981 (see fH)! 



4.1 BBGKY hierarchy 

The Hartree-Fock mean field approximation can be implemented in a variety of ways. We 
have chosen here the approach in terms of the BBGKY hierarchy, truncated to first order. 

4.1.1 Few body-density matrices 

We have already introduced in §^ the concept of the one-body density matrix. We revisit 
here this notion and extend it to two-body density matrices. 

• For a fixed total number of particles 

Let us first consider a system with a fixed total number of particles and let (Xi 2...Ar 
be the -body density matrix. Starting from ai^2...N we introduce simpler objects as 
the one-body and two-body density matrices pi and pu , by taking the trace over the 
states of all the particles but one or two: 

pS^) = iVTr2,3...7v(^i,2,...iv) (98) 
plf = iV(iV -l)Tr3,4...iv(ai,2,...7v). (99) 

In practice the knowledge of pi and pi2 is sufficient to describe most of the experimental 
results. As you know, {r\pi\r) is the density of particles and {fl, r^|pi|rl, r^) is the pair 
distribution function. 
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• For a fluctuating total number of particles 

If N fluctuates according to the probability distribution , we define few-body density 
matrices by the following averages over N : 



Pi = EPnpT^ (100) 

TV 

(N) 



Pl2 = E^ivMf (101) 
N 

Alternatively on can define directly the one-body and two-body density matrices in 
second quantization: 

(rl|pi|r1) = (V^^(ri)V'(r-I)) (102) 
(rl,r^|pi2|r3,rl) = (V''''(^3)V'^(^4)V'(^2)V'(n))- (103) 

Note that the few-body density matrices are normalized as 

Tr[pi] = (N) (104) 
Tr[pi2] = (iV(iV-l)) (105) 

so that one can obtain the variance of the fluctuations in the number of atoms from the 
one-body and two-body density matrices. 

4.1.2 Equations of the hierarchy 

The idea of our derivation of the mean field approximation is to get an approximate closed 
equation for pi by closing the hierarchy with some "cooking recipe" giving pi2 in terms 
of pi . 

To derive the first equation of the hierarchy we start from the exact master equation: 

ih-^cri 2„N = [H, (7i,2..Ar] (106) 
at 

where the Hamiltonian is the sum of one-body and two-body terms: 

N I ^ 

H = Eh^ + ^ E (107) 

i=l i^j-ij=l 



Hartree-Fock approximation 



36 



The single particle Hamiltonian hi contains the kinetic and trapping potential energy of 
the atom i and Vij in the interaction potential between the atoms i and j . Now we 
take the trace of the master equation over the particles 2,3 ... N and multiply it by , 
obtaining 

d f ^ 1 

^h-pl = [/ii,pi] + ATTra,...^ |^^Jl^i,-,(Ti,2....7v] > ■ (108) 

We have kept here only the terms involving the atom 1, as the other terms are commutators 
of vanishing trace. The sum over j amounts to — 1 times the same contribution, e.g. 
the j = 2 contribution, as the atoms are indiscernible. We finally obtain the first equation 
of the hierarchy: 

ihj^pi = [/ii,Pi] + Tr2{[ri2,Pi2]}. (109) 

The equation Eq. ( |109| ) is not closed for pi , as it involves pi2 • The next equation of 
the hierarchy, the equation for pi2 , involves pi23 , etc, up to the -body density matrix, 
where the hierarchy terminates. The mean field approximation consists in replacing pi2 
by an ad hoc function of pi . 



4.2 Hartree-Fock approximation for T > 

4.2.1 Mean field potential for the non-condensed particles 

We use the following simple approximation to break the hierarchy: 

pi2 ^ pf2^ = (^^7^^) Pi ® Pi (^TT^) = + ® 

where P12 is the permutation operator exchanging the states of the particles 1 and 2. 
The last identity in ( |110| ) is obtained by using the commutation of P12 and pi® pi , and 
the fact that = 1 • 

The factorized prescription pi2 = pi® pi is the Hartree approximation. It assumes 
weak correlations between the particles. Indeed at short distances ri2 , the real pi2 
is expected to be a statistical mixture of scattering states of the interaction potential. 
Neglecting the correlations in pi2 between particles 1 and 2 amounts to considering 
only separable, plane wave scattering states, which corresponds to the zeroth order in 
the Born expansion of the scattering theory. Actually pi2 appears in Eq.( [109| ) inside a 
commutator with V12 , so that taking the zeroth order approximation for the scattering 
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states in pi2 corresponds to the first order of tlie Born approximation in tlie equation 
for pi . 

As we are dealing with bosons we have supplemented the Hartree approximation by a 
bosonic symmetrization procedure, involving the permutation operator P12 . Note that 
the symmetrization as it was written works only for particles 1 and 2 in orthogonal states: 



V2 

as the factor is the correct normalization factor only in this case. This is almost true 
for a non-degenerate Bose gas. This restriction forces us to treat separately the case in 
which a condensate is present ( T < ) . 

We now insert the Hartree-Fock ansatz for pi2 in the hierarchy |^ 

zhj^p, = +Tr2{[Vl2,pf/]}. (112) 

In the commutator with Vu we will encounter 

Sin - r^2)(l + P12) = (1 + PuW-l - rl) = 25(rl - r^). (113) 

The fact that P12 commutes with V12 is due to the parity of the delta distribution, 
and P12 acting on a state with two particles at the same position can be replaced by 
the identity. As a consequence, with our zero-range interaction potential, the Fock term 
simply doubles the Hartree term. We finally obtain 

^ (114) 



^ + f/(f) + V(f),pi 
2m 



where V{f) is the mean field potential 

V{f) = 2g{f\pi\f) = 2gp{f). (115) 



The Hartree-Fock Hamiltonian is then 

r2 



h''^il) = ^ + U{f) + 2gpif). (116) 
2m 



■^Note that for the present calculation the regularization of the pseudo-potential is not necessary. 
Indeed by considering plane waves as scattering states in pi2 we suppress any problem of divergences 
in the commutator with V12 , and we can then take V12 as a simple delta distribution. 
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The problem is then formally reduced to the one of an ideal Bose gas moving in a 
self-consistent potential. For g > the mean field corresponds to repulsive interactions, 
as 2gp[f) expels the atoms from the region of high density, while for g < the mean 
field corresponds to attractive interactions. 



4.2.2 Effect of interactions on Tc 

Let us now consider the Hartree-Fock one-body density matrix at thermal equilibrium; 
we use the same formula as the ideal Bose gas Eg. (pQ]) , replacing hi by the Hartree-Fock 
Hamiltonian: 

/5i = {exp[/5(/i^^(l)-M)] -1}-^ (117) 

For ksT ^ AE where where AE is the level spacing of h^^ we can perform the 
semiclassical approximation. We obtain for the spatial density as in Eq. (|55D : 

Psc{r) = -^^73/2(zexp[-/3(f/(f ) + 2(7p,,(f ))]) (118) 

At T = Tc the argument of gs/2 goes to 1 in the point rmin where the potential is 
minimal, so that Einstein's condition still holds in the Hartree-Fock approximation: 

Psc(r„.in)AL = C(3/2). (119) 

For the harmonic trap U{f) = the minimum occurs at the center of the trap, 

^min = so that the chemical potential at the phase transition is given by 

fi = 2gp,ci0). (120) 

It is shifted by the mean field effect with respect to the ideal Bose gas. Using as a small 
parameter p^ci^jg /ksT^ , one can derive at constant the first order change in 

the critical temperature with respect to T° , the transition temperature of the ideal Bose 
gas: 

^ = -2.5pL/^(0)a = -l-33^^-^ArV6. (121) 

For = 10^ atoms of Na in a trap of harmonic frequency u = 27t x 100 Hz, with a 
scattering length a = SOaBohr we find T° ^ Ip K, and 5Tc/T^ — —2.5 x 10^^ , an effect 
for the moment smaller than the experimental accuracy. The fact that 5Tf. is negative 
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for effective repulsive interactions ( a > ) is intuitive: for fixed values of and T the 
interacting gas has a lower density at the center of the trap than the ideal Bose gas, so 
that one needs to further cool the gas to get Bose-Einstein condensation. 

• A calculation of ST^ beyond mean Geld 

The purest situation to study the effect of the interactions on the critical temperature 
is the case of atoms trapped in a flat bottom potential; in this case the density is uniform, 
the previously mentioned intuitive mean field effect is suppressed, and our Hartree-Fock 
theory predicts the same critical temperature as the ideal Bose gas. This prediction 
is actually not correct, and rigorous results for the first order correction of Tc in ap^^^ 
have been obtained recently, by a combination of perturbative theory and Quantum Monte 
Carlo calculations p3| : 

= (2.2 ± 0.25)api/3 + o{ap'/''). (122) 



Recent calculations in the many body Green's function formalism confirm this result ||24 
This effect, if heuristically extended to the trap, is of opposite sign and of the same order 
of magnitude as the mean-field prediction. 



4.3 Hartree-Fock approximation in presence of a condensate 
4.3.1 Improved Hartree-Fock Ansatz 

As already emphasized in the previous subsection the symmetrization procedure of the 



Hartree-Fock prescription Eq. (|110|) has to be modified in presence of a condensate. To 



this end we split the one-body density matrix as 

pi = (Ao)|0)(0|+p; (123) 

where is the condensate wavefunction, (No) is the mean number of particles in the 
condensate and p[ is the one-body density matrix of the non-condensed fraction. The 
Hartree approximation for the two-body density matrix now reads: 



Pi® Pi = (^o)^|0j 01 + remaining Hartree part. 



(124) 
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The first term in the right hand size is aheady symmetrized; the second term can be 
symmetrized as in Eg . ([l 101) as it does not involve coexistence of two atoms in the (only) 
macroscopically populated state (j) . We therefore put forward the following Hartree-Fock 
ansatz: 



— (A^o)^|0) 01 + ( I remaining Hartree part i ~'~ 



Eliminating the remaining Hartree part with the help of Eq.( |124| ), we finally obtain 



(125) 



V V2 J'''-^'\ ^ J v-u/ (126) 

In this way we have avoided the double counting of the condensate contribution that 
would have resulted from the prescription Eq. (|110|) . 



4.3.2 Mean field seen by the condensate 

We replace pi2 in the first equation of the hierarchy by the improved Hartree-Fock ansatz. 
The first bit of the ansatz gives the same result as in the case T > Tc , the second bit 
involves the term: 



Tr2([5(r-I-rl),|0,0)(0, 



n)lM0)(0|]. 



Splitting pi as condensate and non-condensed contribution we arrive at 



^ + Uif) + 2gp{f' 
2m 



{No)g\m\MNom{^\ 



+ 



2m 



+ U{f) + 2gp{f),p[ 



(127) 

(128) 
(129) 



The non-condensed particles still move in the mean field potential 2gp{r) . On the 
contrary the atoms in the condensate see a different mean field potential: 



2gp{r)-g{N,mr)\'' = 2gp'if) + g{Nomf)\' 



(130) 



where p' is the non-condensed density and {No)\(f)\'^ is the condensate density. Q This 

result can be interpreted as follows: An atom in the condensate interacts with non- 

careful reader may argue that we forget here the condition of orthogonahty of the cigenstates of 
p'l to (j) ■ Inclusion of this condition is beyond accuracy of the Hartree-Fock approximation. It will be 
carefully included in the more precise number conserving Bogoliubov approach of 
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condensed particles with the effective couphng constant 2g , and it interacts with another 
particle of the condensate with the effective coupling constant g . 

For repulsive effective interactions { g > 0) this is at the basis of Nozieres'argument 
against fragmentation of the condensate in several orthogonal states: in a box of size L 
in the thermodynamical limit, transferring a finite fraction of condensate particles from 
the plane wave p = to an excited plane wave p = 0{h/L) costs a negligible amount 
of kinetic energy per particle but a finite amount of interaction energy. The transferred 
fraction would indeed be repelled with a stronger amplitude ( 2g rather than g ) by the 
atoms remaining in the condensate. 



4.3.3 At thermal equilibrium 

At thermal equilibrium the one-body density matrix of non-condensed atoms is given 
by the usual Bose distribution for the ideal Bose gas, with the trapping potential being 
supplemented by the mean-field potential: 

Pi = nrrr^ — ...J . n — r (131) 



exp 



2_ + uif) + 2gp{r) - p 



The condensate wave function has to be a steady state of the total, mean field plus 
trapping potential seen by an atom in the condensate: 

M{r) = -^M+[U{f)+ g{No)W)? + 2gp\f)]<j){f). (132) 

The Hartree-Fock single particle energy A should not be confused with the energy per 
particle in the condensate, as it will become clear in the next section. The occupation 
number of the condensate is related to A by the Bose formula: 

We now have to solve in a self consistent way the three equations Eq.( |131|Jl32| , |133|) . 
In practice, when (Aq) is already large, one can assume X = fi , which eliminates one 
unknown A and one equation Eq. (P-33|). 



4.4 Comparison of Hartree-Fock to exact results 
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4.4.1 Quantum Monte Carlo calculations 



N(r)/N 
0.12 

0.10 



0.04 



0.02 



The Quantum Monte Carlo method developed by David Ceperley and others allows to 
sample in an exact way the -body distribution function of a gas of interacting 
bosons at thermal equilibrium. I.e. the algorithm generates random positions fi,...,rN 
for the particles with a probability distribution given by the exact -body distribu- 
tion function of the atoms. 

On the figure ^ the Hartree-Fock prediction for the radial density of particles in a 
spherical harmonic trap, r^p(r) , is compared to the Quantum Monte Carlo result for 
several temperatures below . The Hartree-Fock prediction is in good agreement with 
the exact result, except close to Tc where it tends to underestimate the number of 



particles in the condensate [25 
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Figure 8: Radial density of particles, r'^p{r) , for an interacting Bose gas at thermal equilib- 
rium in an isotropic harmonic trap. Noisy lines: results of a Quantum Monte Carlo simula- 
tion. Smooth solid lines: Hartree-Fock prediction. The curves corresponds to the temperatures 
T /T^ = 0.88 (a), T/T^ = 0.7 (b). The number of particles is = lO'' and the parameters 

87 



are the ones of Rb. These figures are taken from l2 



4.4.2 Experimental results for the energy of the gas 

At JILA the sum of kinetic and interaction energy of the atoms was measured as function 
of temperature, as we have already explained in §^.3|. Whereas the ideal Bose gas model 
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was clearly getting wrong for T < , the Hartree-Fock prediction is consistent with 
the experimental results over the whole considered temperature range (see figure |^). 




Figure 9: Expansion energy of the gas per particle and in units of ksT^ as function of the 
temperature in units of . The filled rhombi correspond to the experimental results of |11]. 



The straight solid line is the prediction of Boltzmann statistics. The dotted curve is the ideal 
Bose gas prediction. The circles are the numerical solution to the Hartree-Fock approach. The 
curved solid line and the dashed line are approximate solutions to the Hartree-Fock equations. 
The inset is a magnification showing the change of slope of the energy as function of T close 
to T = . The figure is taken from |2f:]. 



At very low temperatures ( T < Tc/2 ), measurements at MIT have shown that the 
same energy becomes mainly a function of the number of particles Nq in the condensate. 
By setting p' = in the Hartree-Fock approximation, and using approximations pre- 
sented in the coming section § |5.3| , an analytical expression can be obtained for the energy, 
in excellent agreement with the experimental results (see figure the energy per particle 
has a power law dependence with Nq , with an exponent 2/5 , to be contrasted with the 
constant ideal Bose gas result, and has typical values an order of magnitude larger than 
the zero-point energy of the harmonic oscillator. 
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5 Properties of the condensate wavefunction 



In this section we consider the regime of an almost pure condensate, where the non- 
condensed cloud has a negligible effect on the condensate. At thermal equilibrium with 
temperature T this regime corresponds to the limit T <^ . As we shall see most of the 
experimental results obtained with almost pure condensates can be well reproduced by a 
single equation for the condensate wavefunction, the so-called Gross-Pitaevskii equation, 
derived independently by Gross and Pitaevskii pB . 



5.1 The Gross-Pitaevskii equation 
5.1.1 From Hartree-Fock 

Let us assume that the density of non-condensed particles is much smaller than the density 
of condensate particles over the spatial width of the condensate: 



p'(f)«Aro|0(f)p 



(134) 



where A^^o is the mean number of particles in the condensate, is the condensate wave- 
function normalized to unity: 



(fr(f){f,t)(f)*{f,t) = 1. 



(135) 



In the Hartree-Fock expression of the mean field potential seen by the condensate, derived 
in the previous section we can drop the contribution of the non-condensed particles, 
to get for the evolution of the condensate contribution to the 1-body density matrix: 



th^^{No\<l>){<l>\) 



P 



^ + U{f,t)+gNo\(j){f,t)\\N, 



2m 



(136) 



This equation leads to A'o = constant and to the evolution equation for the condensate 
wavefunction: 



ihdt(p{r, t) 



I^A + f/(f, t) + Nog\<P{f, t) |2 - e(t) 



(f,t). 



(137) 



This non-linear Schrodinger equation is the so-called time dependent Gross-Pitaevskii 
equation. This equation is determined from our Hartree-Fock approach up to an arbitrary 
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real function of time, ^{t) , as Eq.( p.36| ) involves a commutator to which ^(t) does not 
contribute. In general the precise value of ^{t) is considered as a matter of convenience, 
as it can be absorbed in a redefinition of the global phase of . The knowledge of 
the value of ^(t) can become important when one is interested in the evolution of the 
relative phase of two Bose- Einstein condensates. The value of ^(t) has been derived in 



29[| assuming a well defined number of particles in the condensate. If the condensate is 



assumed to be in a Glauber coherent state that is a quasi-classical state of the atomic 
field with a well defined relative phase (see §|D one obtains C,{t) = as we will see in 



When the gas is at thermal equilibrium, the only time dependence left for is a 
global phase dependence. The most convenient choice is to assume dt4> = so that ^(t) 
is a constant. As shown in § 4.3.3 this constant is very close to the chemical potential 
of the gas as Nq is large so that we get the so-called time independent Gross-Pitaevskii 
equation: 



-—A + U{f) + Nog\(j){f)\' 
Zm 



(138) 



Both the time independent and the time dependent Gross-Pitaevskii equations can be 
solved numerically. But, as explained in the next part of this section, the fact that the 
trap is harmonic allows one to find very good approximate analytical solutions. 



5.1.2 Variational formulation 

Variational calculus turns out to be a very fruitful approximate technique in the solution 
of the Gross-Pitaevskii equation. We therefore derive here a variational formulation of 
the Gross-Pitaevskii equation. 

• Time independent case 

The time independent Gross-Pitaevskii equation can be obtained from extremalization 
over of the so-called Gross-Pitaevskii energy functional: 

E[<f),<p*] = No Jd^f 



h 



1 



— igrad0|2 + u{fmn\' + ^Nogmni^ 



(139) 



with the constraint that is normalized to unity. 
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Proof: We take into account the normalization constraint with the method of La- 
grange multipher, so that we simply have to express the fact that (p extremalizes without 
constraint the functional: 



X[(l>,f] = E[(f>,(f>*] - XNo I (fr(t>{r)(t>*{r). 



(140) 



The parameter A is the Lagrange multiplier. We calculate the first order variation of X 
due to an infinitesimal arbitrary variation of the condensate wavefunction: 



We obtain: 



2m 



(f){r) ^ (f){f) + S(f){f). (141) 



grad 50* • grad + U{r)5(t)*(t) + Nog5(l)*(f)*(f)'^ - \5(t)*(t) + c.c. (142) 



We modify the variation of the kinetic energy term by integrating by part, assuming that 
vanishes at infinity: 



j d^f (grad Sep* ■ grad + c.c.) ^ - j d^f {S(j)*A(f) + c.c). 



(143) 



The variation SX has to vanish for any 6(j) . We can take as independent variables the 
real part and the imaginary part of 6(f) , or equivalently 6(f) and 6(f)* as it amounts 
to considering independent linear superpositions of the real and imaginary part. We 
therefore obtain: 



2m 



A + U{f) + Nog\(l){f)\'-X 



0. 



(144) 



We recover the time independent Gross-Pitaevskii equation, with X — /i , which gives a 
physical interpretation to the Lagrange multiplier A . 

• Time dependent case 

The time dependent Gross-Pitaevskii equation with the choice ^{t) = is obtained over 
a time interval [ti, ^2] from extremalization of the action: 



dt 



|[(0|||0)-c.c.]iVo-£;[0(i),0*(t)] 



(145) 



with fixed values of |0(t = ti)) and \(f){t = ^2)) • 
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• Physical interpretation of the Gross-Pitaevskii energy functional 

We now show that E[(f),(j)*] is simply the mean energy of the gas in the Hartree-Fock 
approximation in the hmit of a pure condensate. As the -body Hamiltonian is a sum 
of one-body and two-body (binary interaction) terms, 

TV ^ 

H = T.h^ + ^T.V^J (146) 

1=1 iytj 

the mean energy of the gas involves the one-body and two-body density matrices: 

{H) = Tr[/iipi] + ^Tr[Vi2Pi2]. (147) 
In the limit of a pure condensate we keep only the condensate contribution to pi : 

pi~iVo|</')(</'| (148) 
and we approximate pi2 by the Hartree ansatz 

pi2~piOpi. (149) 

We then obtain E[(f),(j)*] = (H) . It was actually clear from the start that E[(f),(j)*] was 
the sum of kinetic energy, trapping potential energy and mean field interaction energy of 
the condensate. 

A different and interesting point of view at zero temperature is to use directly a 
Hartree-Fock ansatz for the ground state wavefunction |\E') of the gas, assuming that all 
the particles are in the condensate: 

|$) = |iV:0) = |0)®...®|0). (150) 

The mean energy of |\E') for the interaction potential g6{fi — r2) is then 

= N j d?r 

which differs from Eq. (|139| ) in the limit Nq = N only by the occurrence of a factor (A^— 1) 
rather than A^ in front of the coupling constant g , ensuring that the interaction term 
disappears for A^ = 1 ! 



2m 



|grad0|2 + f/(f )|0(f )|2 + -{N - l)g\(p{f)\' 



(151) 
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• What is the chemical potential ? 

At zero temperature, assuming a pure condensate Nq c:^ N , the usual thermodynamical 
definition of tlie cliemical potential fi reduces to: 

where we have made appear the explicit dependence of E on Nq . When one takes the 
total derivative of E with respect to A^^o , one gets in principle a contribution from the 
implicit dependence of E on Nq through the Nq dependence of 0, 0* ; actually this 
contribution vanishes as the variation of E due to a change in 0,0* vanishes to first 
order in this change. We therefore get 

|-|gr;d,#.p + U{rMf)\' + N„g\4,{f)\' 



I 



2m 



(153) 



This quantity coincides with the chemical potential indeed, as can be checked by multi- 
plying the time independent Gross-Pitaevskii equation by 0* and integrating over the 
whole space. As g does not have the factor 1/2 in Eq.( |153| ), whereas it is multiplied by 
1/2 in the expression for -E[0, 0*] , we see that in the interacting case g ^ Q : 

, E , , 

^ TT (154) 
iVo 

that is the chemical potential /i differs from the mean energy per particle. 

5.1.3 The fastest trick to recover the Gross-Pitaevskii equation 

Starting from the second quantized form of the Hamiltonian, 

H = [ dl^\l)hMl) + l [ dlf rf2^t(i)^t(2)v^^2^(2)^(l) (155) 



2 

where 1 and 2 stand for three-dimensional coordinates in real space, one first derives 
the Heisenberg equation of motion for the field operator: 

^hJ^^Pil) = [i,{l),H] = d^,^,^H (156) 
= h^i}{l)+ j d2i,\2)Vi2i^{2)ij{l) (157) 
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and then replaces the quantum field operator by a classical field: 




As Vi2 is the pseudo-potential, the equation that we get for (f) is the time dependent 
Gross-Pitaevskii equation with ^(t) = . 

This sheds a new light on the Gross-Pitaevskii equation: the Gross-Pitaevskii equation 
is the equation of motion of the atomic field in the classical approximation, neglecting 
quantum fiuctuations of the field. A Bose-Einstein condensate is a classical state of the 
atomic field, in a way similar to the laser being a classical state of the electromagnetic 
field. 



5.2 Gaussian Ansatz 

In this subsection we look for a variational solution to the Gross-Pitaevskii equation in a 



harmonic trap, using a Gaussian ansatz for [^. The choice of a Gaussian is natural 
in the non-interacting limit g ^ , where it becomes exact. It turns out to give also 
interesting results in presence of strong interactions. 



5.2.1 Time independent case 

Consider for simplicity an isotropic harmonic trap, where the atoms have the oscillation 
frequency u . We assume the following Gaussian for the condensate wavefunction: 



1 



the spatial width a being the only variational parameter. The mean energy per particle 
can be calculated exactly for this ansatz: 

E[<P,<P*] 3 2 , h^Noa 1 

e = — — = + -mw a H — 1=. (160) 

A^o 4ma2 4 m V27r 



The form of the result is intuitive: the kinetic energy term scales as Ap^. , where Ap^ = 
h/{2Ax) = h/ (\/2(t) ; the trapping potential energy scales as o"^ and the interaction en- 
ergy per particle is proportional to the coupling constant g = Anfi^a/m and to the typical 
density of atoms in the gas, N^/a^ . Taking the harmonic oscillator length (h/mcuY^'^ 



The condensate wavefunction 



50 



as a unit of length and the harmonic quantum of vibration huj as a unit of energy we 
get the simple form: 

^ (161) 



3 



1 2 



2a3 



where the only physical parameter left is 

This parameter x nieasures the effect of the interactions on the condensate density: The 
case X ^ 1 corresponds to the weakly interacting regime, close to the ideal Bose gas 
limit X = ; the case x ^ 1 corresponds to the strongly interacting regime. 

• case a > 

In the case of effective repulsive interactions between the particles, the dependence of e 
with a is plotted in figure [1^. In the limit a — , the energy e is dominated by the 
positively diverging repulsive interaction ( ~ 1/cr^ ). For large a the trapping potential 
term ~ cr^ dominates. The function e has a single minimum, in a = ctq , solving 

de 

— {a^) = Q^al = a^ + X- (163) 

For X ^ 1 recovers the ground state of the harmonic trap, with ao = 1 . For x ^ 1 
the condensate cloud becomes much broader than the ground state of the harmonic trap, 

ao - x}'^ oc nI'^. (164) 

In this regime one can check that the kinetic energy term becomes negligible as compared 
to the trapping energy: 

-Enn 11 . . 

= (165) 



SO that the steady state of the condensate is an equilibrium between the trapping potential 
and the repulsive interactions between particles. This regime will be studied in detail in 
the next subsection. 



• case a < 
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Figure 10: Energy per particle in the condensate in units of Hlo as function of the variational 



For effective attractive interactions between the particles the shape of e as function of 
a depends on the balance between kinetic and interaction energy (see figure |lTD. The 
interaction energy is negatively diverging as a ^ always faster than the positively 
diverging kinetic energy so that cr = is always a minimum of e , with e = —oo : 
the condensate is in a spatially collapsed state ! Of course the Gross-Pitaevskii equation 
not longer applies for a too small a , as the validity of the Born approximation requires 
fc|a|~|a|/o"<^l. For |x| larger than some critical value \Xc\ , this collapsed minimum 
is the only one of e{a) so that we do not find any stable solution for the condensate 
wavefunction. When |x| is smaller than \xc\ the kinetic energy term, which is opposed 
to spatial compression of the gas, is able to beat the attractive energy over some range of 
a , so that a local minimum of e{a) appears, in a = ao , separated from the collapsed 
minimum by a barrier. 

To calculate \xc\ we express the fact that the stationary point of e in a = ao has 
now a vanishing curvature (infiexion point of £ ): 



width a in units of (h/muj)^^'^ . Case of effective repulsive interactions a > . 








(166) 
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Figure 11: For an effective attractive interaction a < between the particles energy per 
particle in the condensate in units of huj as function of the variational width a in units of 
(n/mL<j)i/2 . The curve has two possible shapes, (a) with two minima when |x| is smaller than 
a critical value |xc| > and (b) with a single minimum for \x\ > \Xc\ ■ 



X=Xc 







(167) 



By eliminating ao between these two equations we obtain 

4 



Xc 



55/4 



-0.5350... 



;i68) 



This result can be rephrased in terms of a maximal number of atoms A^^q that can be 
put in the condensate without inducing a collapse, according to a Gaussian ansatz: 



NS\a\ 



fi/muj 



\Xc 



-0.67. 



(169) 



A more precise result has been obtained by a numerical solution of Gross-Pitaevskii 
equation, not restricting to the subspace of Gaussian wavefunctions no solution of 
the time independent Gross-Pitaevskii equation is obtained for Nq > Nl^ , where 



N^\a\ 
h/muj 



-0.57. 



(170) 
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By a generalization of the Gaussian ansatz to the case of a non-isotropic harmonic 
trap one can also get a prediction of Nq for the parameters of the lithium experiment 



of Hulet's group [32]. In the experiment the traps frequencies are Uz = 2n x 117 Hz and 
^x,y = 27T X 163 Hz, and the scattering length is a = —27 Bohr radii. The Gaussian 
prediction is then ~ 1500 , consistent with the experimental results. 

• Physical origin of the stabihzation for a < 

In a harmonic trap, the energy of the ground state level is separated from the energy 
of excited states by hu . At low values of x the mean interaction energy per particle, 
~ p\g\ , where p is the density, is much smaller than hu so that it cannot efficiently 
induce a transition from the ground harmonic level to excited harmonic levels. Initiation 
of collapse on the contrary requires that the wavefunction (p expands on many excited 
levels in the trap, so that the density |0p can exhibit a high density peak narrower than 
h/muj . We therefore intuitively reformulate the non-collapse condition as 

p\g\ < huj. (171) 



Estimating p as Nq/ {Ti/muj)'^/'^ we recover a A^q scaling as \Jh/muj/\a\ . This reasoning 
also applies to the gas confined in a cubic box with periodic boundary conditions, as we 
shall see in section §H of the lecture. 

5.2.2 Time dependent case 



As done in [^, ^ the Gaussian ansatz can be generalized to the time dependent case. 
We assume here for simplicity that the condensate, initially in steady state, is excited only 
by a temporal variation of the trap frequencies uJa{t) ; then no oscillation of the center 
of mass motion of the condensate takes place, remaining of vanishing mean position 
and momentum. The Gaussian ansatz then contains only exponential of terms quadratic 
with position, its does not involve exponential of terms linear with position: 



(172) 



We do not assume that the trap is isotropic, so we have as variational parameters 3 spatial 
widths cTq, { a = x,y, z ), 3 factors 7q, governing the spatially quadratic phase and a 
global phase 6 . 
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One gets time evolution equations for the variational parameters by inserting the 
ansatz for in the action A of Eq. (|145D and by writing the Lagrange equations ex- 
pressing the stationarity condition. It turns out that ■ja can be expressed in terms of 
the widths and their time derivatives: 

so that one is left with equations for the a a 's. Taking u!~^ as a unit of time, ^Jh/mu 
as a unit of length, where u is an arbitrary reference frequency, we get: 

<Ta + l^l(Ta = ^ + (174) 

0"a O'aO'xO'yO-z 



where the trap frequencies are Ua = Va'^ and x is defined in Eq.( [162|) . In the absence 
of interaction ( x = ) these evolution equations become exact, and give a remarkable 
(and known !) result for the time dependent harmonic oscillator. In the interacting case 
( X 7^ ) these equations can be cast in Hamiltonian form as the "force" seen by the 
variable cTq derives from a potential. The corresponding dynamics is non linear and non 
trivial; chaotic behavior has been obtained in |^ in the limiting regime of x ^ 1 where 
the l/o"^ can be neglected. 

One can use Eq. (|174|) to study the response of the condensate to a weak excitation. 



the trap frequency lo^ in the experiments being typically slightly perturbed from its 
steady state value Wq,(0) for a finite excitation time. Linearizing the evolution equations 
in terms of the deviations of the a 's from their steady state value: 

(y^{t) = at^ba^{t) (175) 

one gets a three by three system of second order differential equations for the 6a 's. 
Looking for eigenmodes of this system, one finds three eigenfrequencies [0]. Their values 



have been compared to experimental results at JILA ||36[, see Fig.|T2|: the agreement is 
very good, not only in the weakly interacting regime x ^ 1 but also in the regime 
X ^ 1 , where the Gaussian ansatz for the condensate wavefunction has no reason to be 



a good one! The explanation of this mystery is given in §[5.4.1 
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Figure 12: Frequencies of two eigenmodes of a condensate in a cylindrically symmetric harmonic 
trap, in units of the radial trap frequency , as a function of a parameter proportional to x 



measuring the strength of the interactions. Plotting symbols: measurements at JILA |36|. Solid 



lines : predictions of the Gaussian Ansatz |34] 



5.3 Strongly interacting regime: Thomas-Fermi approximation 

In this subsection we focus on the strongly interacting regime: the scattering length is 
positive, with the dimensionless parameter x of Eg. ( |162| ) much larger than one. This 
regime is the so-called Thomas-Fermi regime. As we now see analytical results can be 
obtained in this limit. 



5.3.1 Time independent case 

If we put a large enough number of particles into the condensate the atoms will experience 
repulsive interactions that will increase the spatial radius of the condensate to a value R 
much larger than the one of the ground state of the harmonic trap: 

/ h 

> . (176) 
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For increasing value of A^^q , R increases so that the momentum width of the condensate, 
scahng as h/R as 0o is a non-oscillating function of the position, is getting smaller and 
smaller. More precisely we find that the typical kinetic energy of the condensate becomes 
much smaller than the typical harmonic potential energy of the condensate: 

^ ^ — ^ <1. (177) 

-E'harm muo'^R^ \muR'^j 

The mechanical equilibrium of the condensate in the trap then comes mainly from the 
balance between the expelling effect of the repulsive interactions and the confining effect 
of the trap. 

In this large R regime we neglect the kinetic energy term in the Gross-Pitaevskii en- 
ergy functional, which leads to a functional of the condensate density only (similarly to the 
Thomas- Fermi approximation for electrons). This approximation amounts to neglecting 
the term in the Gross-Pitaevskii equation: 

/i0(f) ~ f/(f)0(f ) + A/'o^|0(f )|V(r ). (178) 

Taking to be real we find that 



fi — U{r)\ 



(179) 



V ^0^7 ; 

in the points of space where /i > U{f) , otherwise we have 0(r ) = . 

This very important, Thomas- Fermi result Eq.( |179| ) can also be obtained in a local 
density approximation point of view. A spatially uniform condensate with a chemical 
potential fi and in presence of a uniform external potential U has a density A'^ol^P = 
{fi — U)/g. Applying this formula with a r dependent potential U gives again Eq. (|T7Up . 
A local density approximation can be used only if the density of the condensate varies 
slowly at the scale of the so-called "healing length" ^ , introduced in § ^.3.4| ; one can check 
that the condition ^ <^ i? is indeed satisfied in the Thomas-Fermi regime. 

We specialize Eq. ( |1 79|) to the case of a harmonic but not necessarily isotropic trap: 

U{f) = l:mY^u;lrl (180) 

where a = x,y, z label the eigenaxis of the trap. The boundary of the condensate 
fi = U{f) is then an ellipsoid with a radius Ra along axis a given by: 

/X = ImcolRl (181) 
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The condensate wavefunction can be rewritten in terms of these radii: 



Nog. 



1/2 




(182) 



Using the normalization condition of (f) to unity we can also express the "normalization" 
factor yJ^i/Nog in terms of the radii. The integral of |0p can be calculated in spherical 
coordinates after having made the change of variable Ua = fa/Ra ■ This leads to 



1/2 



/ 



15 



1/2 



V 



571 



;i83) 



Eliminating in terms of /i thanks to Eq. (|181|) we can calculate the chemical potential: 

2/5 



-hid 



15- 



{h/muj 



-\l/2 



where u is the geometrical mean of the trap frequencies: 

a/3 



We can now see that in the limit x ^ 1 the chemical potential /i satisfies 



;i84) 



(185) 



(186) 



which is a convenient way of defining the Thomas- Fermi regime. 

We can now compare these Thomas- Fermi predictions to the MIT experimental results 
on the energy of the condensate [l^. In the experiment the trapping potential is switched 
off abruptly, so that the energy of the gas abruptly reduces to i^red = -^kin + -^int — -^int ; 
afterwards the cloud ballistically expands, Ei^t is converted in kinetic expansion energy 
that can be measured. In the Thomas- Fermi approximation the integral of NQg\(l)\'^/2 
can be done, which leads to 



E\r 



7 



iVo/i oc nI'\ 



(187) 



The resulting dependence in A^^o is in good agreement with the MIT results, see Fig.|^. 

From the expression of the chemical potential we can also calculate the total energy 
of the condensate in the trap , as /i = OnqE : integrating over A^o gives 



E ^ -fxNo. 

One can then check explicitly that fi ^ E/NqI 



(188) 
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5.3.2 How to extend the Thomas-Fermi approximation to the time dependent 
case ? 

We would like to analyze time dependent situations encountered in the experiments, e.g. 

• ballistic expansion of the gas: this is a crucial example, as it is a standard experi- 
mental imaging technique of the condensate 

• collective excitations: response of the condensate to a modulation of the trap fre- 
quencies 

in the strongly interacting regime. An immediate generalization of the Thomas-Fermi 
approximation consisting in neglecting the kinetic energy of the condensate is now too 
naive! In the case of ballistic expansion for example the interaction energy is gradu- 
ally transformed into kinetic energy when the cloud expands so kinetic energy becomes 
important! 

The trick is actually to split the kinetic energy in two contributions, one of them 
remaining small and negligible in the time dependent case. This is performed using the 
so-called hydrodynamic representation of the condensate classical field, split in a modulus 
and a phase: 

N^/^(l){r) = p^/\r)e'^^'^/^ (189) 

where S has the dimension of an action and p is simply the condensate density. The 
mean kinetic energy of the condensate then writes 



^kin[0,0*] = J — Igrad 



x2 (grad^) 
— grad^/p +P 



2m ^ ' 2m 



(190) 



As we shall see during ballistic expansion of the condensate the density p remains a 
smooth, slowly varying function of the position so that it has a very small contribution to 
the kinetic energy; most of the kinetic energy induced from interaction energy is stored 
in the spatial variation of the phase of the condensate wavefunction. 
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5.3.3 Hydrodynamic equations 



In this subsection we rewrite the time dependent Gross-Pitaevskii equation in terms of 
the density p and the phase S . This can be done of course by a direct insertion of 



Eq.( [L89| ) in the Gross-Pitaevskii equation. 

A more elegant way is to use the covariant nature of the Lagrangian formulation of 
the Gross-Pitaevskii equation, Eq.( 145 ). We rewrite the density of Lagrangian in terms 
of p and S : 

/ ^ \2 (gradS*) Q „ 

pdtS + — {&^d^) +p^^^—L + U{T,t)p + y . (191) 



An evolution equation for an arbitrary coordinate Q{f,t) of the field is obtained from 
the Lagrange equation: 

dC \ f dC \ d£ 



dt 



(192) 



We first specialize the Lagrange equations to the choice Q = ^ ; dividing the result- 
ing equation by 2y/p we obtain 

(193) 



2m 

Then we set Q = S in the Lagrange equations which leads to 

P 



dtp + div 



m 



-grad S 



0. 



(194) 



This last equation looks like a continuity equation. This is confirmed by the following 
physical interpretation of grad S . It is known in basic quantum mechanics that the 
probability current density associated to a single particle wavefunction is 

h 



Jproba 



2im 



*grad ( 



c.c. 



(195) 



Multiplying this expression by A^q ? as there are Nq particles in the condensate, and 
introducing the (p, S) representation of we get the following expression for the current 
density of condensate particles: 



J = P- 



grad S 



m 



pv 



(196) 
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where v is the so-called local velocity field in the gas. 

Equation ( p,94| ) is therefore the usual continuity equation: 

dtp + div[pv]=0. (197) 

The other equation ( |193| ) can be turned into an evolution equation for the velocity field 
by taking its spatial gradient: 



mdtv + grad 



-mv + U{f) + gp{f) 2^ ^ ^ 



0. (198) 



This looks like the Navier-Stockes equation used in classical hydrodynamics, in the 
limiting case of a fiuid with no viscosity. The term grad(|mf ^) looks unusual but using 
the fact that v is the gradient of a function S/ m one can put it in the usual form of a 
convective term: 

grad ( -mt>^ ) = m{v ■ grad)-!/. (199) 



,2 

A difference with classical hydrodynamics is the so-called quantum pressure term, involv- 
ing 

- 7^ 200 

the only term in the equations ( p.97|jl98| ) where h appears. 



5.3.4 Classical hydrodynamic approximation 

The classical hydrodynamic approximation consists precisely in neglecting the quantum 
pressure term Eq. (|200|) in the equation ( [L98| ) for the velocity field of the condensate. 



We can estimate simply the validity condition of this approximation. Denoting d 
a typical length scale for the variation of the condensate density p(r ) we obtain the 
estimate 

Comparing the quantum pressure term Eq. (|200|) to the classical mean field term pg yields 
the condition 

< 9p{r) ~ Pmaxfl' (202) 



md? 
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where Pmax is the maximal density (usually at the center of the trap). This validity 
condition can be reformulated in terms of the healing length, 

-^^ ■ (203) 

Note that ^ is sometimes also called coherence length, which can be confusing. 

Why this name of healing length for ^ ? Imagine that you cut with an infinite wall a 
condensate in an otherwise uniform potential. Right at the wall the condensate density 
vanishes; far away from the wall the density of the condensate is uniform. The condensate 
density adapts from zero to its constant bulk value over a length typically on the order 
of ^ . This can be checked by an explicit solution of the Gross-Pitaevskii equation: 

N'fy{x, y, z) = p]/l tanh (^-|- j (204) 

where 2; = is the plane of the infinite wall. This explicit solution shows that at 
a distance z ^ ^ from the infinite wall there is no more any effect of the boundary 
condition 0(x, y.z = {)) = . This is to be contrasted with the case of the ideal Bose 
gas: the ground state between infinite walls separated by the length L then scales as 
sm.{Tiz/L) and depends dramatically on L. 

For a moderate excitation of the condensate by a modulation of the trap frequencies, 
or in the course of ballistic expansion of the condensate, we shall see that the only typical 
length scale for the variation of the condensate density is the radius R of the condensate 
itself. One can then check that in the Thomas-Fermi regime the classical hydrodynamic 
approximation indeed applies: 

^ ^ X ' =T^» 1- (205) 



In the Thomas-Fermi regime we therefore neglect the quantum pressure term to obtain 
dt+{v- grad)] v{f, t) = -grad ([/(f, t) + gp{f, t)) = F{f, t). (206) 



m 



This equation is then a purely classical equation, Newton's equation in presence of the 

— * 

force field F written in Euler's point of view. The operator between square brackets is 
simply the so-called convective derivative. 
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It is instructive to rewrite Eq. (p06| ) in Lagrange's point of view. One then follows a 
small piece of the fluid in course of its motion. Denoting r{t) the trajectory of the small 
piece of fluid we directly write Newton's equation: 

m^f (t) = F{f{t) ,t) = - [grad {U + gp)\ (f (t) , t) . (207) 

This equation automatically implies the continuity equation ( |197| ) and the Euler equation 
( |206|) . The unusual feature is that the force field depends itself on the density of the 



gas, so that we are facing here a self-consistent classical problem, corresponding formally 
to the mean field approximation for a coUisionless classical gas! A surprising conclusion, 
knowing that we are actually studying the motion of a Bose-Einstein condensate! 

5.4 Recovering time dependent experimental results 
5.4.1 The scaling solution 



It turns out that the self-consistent classical problem Eg. ( |207|) can be solved exactly for 



the particular conditions of a gas initially at rest and in a harmonic trap. 

At time t = we assume a steady state Bose-Einstein condensate in the trap, of 
course in the Thomas-Fermi regime so that the classical hydrodynamic approximation 
is reasonable. The steady state of Eq.( |207D corresponds to a force field F vanishing 
everywhere, so that 

U{f) + gp{f) = constant. (208) 

One recovers the stationary Thomas-Fermi density profile, the constant being determined 
from the normalization condition of p and therefore coinciding with the Thomas-Fermi 
approximation for p . 



At time t > the trapping potential remains harmonic with the same eigenaxis 
but the eigenfrequencies of the trap can have any time dependence: 

Uir,t) = ^ Y: mulityl (209) 

a=x,y,z 

Then any small piece of the fluid with initial positions ra{0) along axis a will move 
according to the trajectory 

rait) = A„(t)r„(0) (210) 
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where the scahng factors Xa{t) depend only on time, not on the initial position of the 
small piece of fluid. In other words the density of the gas will experience a mere (possibly 
anisotropic) dilatation 



We can see simply why the ansatz Eq. (|210|) solves indeed Eq. (|207|) for a harmonic trap. 



As the initial density in the trap has a quadratic dependence on position, so will have the 
density at time t . The gradient —gTad{pg) appearing in the expression of the force field 
will then be a linear function of the coordinates; so is the harmonic force — grad ?7(r, t) . 
Newton's equation is therefore linear in the coordinates; dividing it by r^iO) one then 
gets equations for Xait) irrespective of the initial coordinates ra(0) ! 

More details are given in [^8|, ^ , we give here the equations for the scaling parameters: 

i-2Ut)= f -ujl{t)X^{t), a = x,y,z (212) 

with the initial conditions 

XaiO) = (213) 

|a.(0) = (214) 

since the condensate is initially at rest. 

Finally we make the connection between these scaling solutions and the equations for 
the spatial widths cTq obtained in Eq. (|174| ) from a time dependent variational Gaussian 



ansatz for the condensate wavefunction. We are here in the Thomas- Fermi regime x ^ 1 
so that the l/cr^ terms can be neglected in Eq.( |174| ). The steady state solutions for the 
cTq 's are then cTq(O) ~ x^/^z/^/^/z/a(0) where u is the geometrical mean of the initial 
frequencies z/a(0) , and the quantities aa(t)/aa{0) then obey the same equations as the 
Xa 's! The Gaussian ansatz, which has the wrong shape in the Thomas-Fermi regime, is 
however able to capture the right scaling nature of the solution! This explains why the 
collective mode frequencies obtained from Eq. ( |1 74| ) are a good approximation, not only 
when X ^ 1 , where the Gaussian ansatz was expected to hold, but also in the strongly 
interacting regime x ^ 1 • 
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5.4.2 Ballistic expansion of the condensate 

At time t = 0^ the trapping potential is turned off suddenly. The scaling parameters 
then satisfy the simpler equations 

These equations are still difficult to solve analytically. In the experimentally relevant 
regime of cigar-shaped traps, with lJz{0) <^ uj^iO) = ujy{0) , one can find an approximate 
solution H. 

Experimentally the scaling predictions have been tested carefully. First one can see if 
the ballistically expanded condensate density has indeed the shape of an inverted parabola 
||38|| . Second one can measure the radii of the condensate as function of time to see if 
they fit the scaling predictions [0. Both tests confirm the scaling predictions in the 
Thomas- Fermi regime. 

5.4.3 Breathing frequencies of the condensate 

A typical excitation sequence of breathing modes of the condensate proceeds as follows. 
One starts with a steady state condensate in the trapping potential. Then one modulates 
one of the trap frequencies for some finite time texc , at a frequency close to an expected 
resonance of the condensate. Then one lets the excited condensate evolve in the unper- 
turbed trap for some adjustable time tosc • Finally one can perform imaging of the cloud, 
e.g. by performing a ballistic expansion of the condensate and measuring the aspect ratio 
of the expanded cloud. By repeating the whole sequence for different values of tosc one 
can reconstruct the aspect ratio as function of tosc • 

Such a procedure has been used at JILA and at MIT. In figure |1^ are shown results 
obtained at MIT in a cigar-shaped trap, for a modulation of the trap frequency along 
the slow (that is weakly confining) axis z . The solid line corresponds to the prediction 
of scaling theory, the input parameters being (i) the oscillation frequencies Ua of the 
atoms in the trap, (ii) the precise way the excitation is performed, and (iii) the duration 
of ballistic expansion; the agreement between theory and experiment is good, considering 
the fact that there is no fitting parameter . 



If one is interested only in the frequencies of the breathing eigenmodes of the conden- 
sate it is sufficient to linearize the equations of the scaling parameters around their steady 
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state value: 



/3 



(216) 



and find the eigenvalues of the corresponding three by three linear system (see also § |6.3.3|) . 
For a trap with cylindrical symmetry one gets the eigenfrequencies = \/2uj±{0) and 



3cj^(0) + AujKo) ± (9cj^(0) + 16cji(0) - 16c<;^(0)cj1(0) 



1/2- 



(217) 



The mode observed at MIT corresponds to the — sign in the above expression; as the 
trap was cigar-shaped in the experiment, ^^(O) ^ ^^^(O) so that one has the approximate 
formula 

n ^ (^-j ^,(0). (218) 
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Figure 13: Aspect ratio of the excited and ballistically expanded condensate as a function of 
free oscillation time tosc • The expansion time is 40 ms, the unperturbed trap frequencies are 
"^±(0) = 2iT X 250 Hz, LOz{0) = 27r x 19 Hz. Solid line: theory. Diamonds: experimental data 
obtained at MIT. 
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6 What we learn from a linearization of the Gross- 
Pitaevskii equation 

There are several important motivations to perform a linearization of the Gross-Pitaevskii 
equation around a steady state solution 0o : 

• as the Gross-Pitaevskii equation is a non-linear equation it is crucial to check the so- 
called "dynamical" stability of the steady state solution. More precisely one has to 
check with a linear stability analysis that any small deviation 6(f) of the condensate 
wavefunction from (pQ does not diverge exponentially with time. Otherwise 0o 
may not be physically considered as a steady state as even very small perturbations 
will eventually induce an evolution of the condensate wavefunction far from 0o • 

• as a byproduct of linear stability analysis we obtain a linear response theory for the 
condensate very useful to interpret experiments which apply a weak perturbation 
to the condensate. 

• another important byproduct is the Bogoliubov approach which gives a description 
of the state of the non-condensed particles that is still approximate but more ac- 
curate at low temperature (typically < fi ) that the Hartree-Fock approach. 
This allows to check the so-called "thermodynamical stability" of the condensate 
and will be the subject of 

6.1 Linear response theory for the condensate wavefunction 

6.1.1 Linearize the Gross-Pitaevskii solution around a steady state solution 

Let (poir) be a steady state solution of the Gross-Pitaevskii solution in the time inde- 
pendent trapping potential Uo{f) : 



The trapping potential is then slightly modified by a time dependent perturbation 
5U{f, t) , resulting in a total trapping potential 




(219) 



U{f,t) = Uo{f) + SU{f,t). 



(220) 
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The condensate wavefunction, initially equal to 0o ; evolve according to 



ihdf 



- A + U + gNo\<j)\'-fi 
2m 



(221) 



As SU is so small we assume that experiences only a small deviation from 0o • 

(j)if,t) = Mr) + Hir,t) (222) 
. Neglecting the second order product of 



so that we can linearize Eq.( p21| ) in terms of 
66 and 6U we obtain: 



ihdt 



+ 2gNocP*ocl>oS^ + gNo^l5(P* + SUcPq. (223) 



Note the presence of the factor 2 in front of the term proportional to g6(f) ; it turns out 
(and this should become clear in the Bogoliubov approach) that this factor 2 has the same 
origin as the one in the Hartree-Fock potential Eg . (11151) for the non-condensed particles. 
As remains normalized to unity, as 0o "was, we note that to first order in 50 , 



£r [50(f,t)0*(f ) + 0o(f )(50*(f,t)] = 0. 



(224) 



A peculiar feature of Eq. ( f^23| ) is that, though it is obtained from a linearization pro- 
cedure, it is not a linear equation for 50 in the strict mathematical sense: if 50 is a 
particular solution of the homogeneous part of this equation (set 6U = ), the function 
a6(j) (where a is a constant complex number) is generally not a solution of the homoge- 
neous part anymore because of the coupling of 50 to 50* . There are several possibilities 
to restore this linearity. A first one is to consider as unknown functions the real part 
and the imaginary part of 50 . A second, more elegant method, more common in the 
literature, is to introduce formally as unknown the two-component column vector: 



50(f,t) 
50*(f,t) 



(225) 



the functions 50 and 50* being now considered as independent. We then rewrite 
Eq.( p23| ) as the linear system: 



ihdf 



6<j){f,t) 
6(f)* {f,t) 



GP 



50(f,t) 
50*(f,t) 



+ 



S{f,t) 
-S*{f,t) 



(226) 
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with a source term S{f,t) = 6U{f,t)(f)Q{f) and a linear operator 



C 



Hop + gNo\(Po\^ gN^^^ 



(227) 



-gNo(l)f -[HGP + gNo\(j)o 
where we have introduced the Gross-Pitaevskii Hamiltonian: 

Hop = -—A + Uo + gNo\(Po\^ - ft. (228) 
2m 

Note the presence of complex conjugation in the second line of Cgp ', it also applies to 
the potential Uq , without effect here as Uq , hermitian function of r , is real; it should 
not be forgotten if situations where the potential contains a complex term such as —QL^ 
where is the angular momentum operator (inertial term in a frame rotating at angular 
velocity Q ). 

As the operator Cgp is time independent the general method to determine the time 
evolution of 6(j) is to diagonalize Cgp and expand 6(f) on the corresponding eigenmodes. 
At this stage one faces a slight difficulty: it turns out that Cgp is not diagonalizable, 
that is the set of all eigenvectors of Cgp does not form a basis (in general one vector 
is missing to span the whole functional space). Mathematically this can be solved by 
putting Cgp into the so-called Jordan normal form. Here we use the more physical 
following procedure. 

6.1.2 Extracting the "relevant part" from S(f) 

We split (50 on a component along 0o and a part orthogonal to 0o • 

Scj>{f,t) = vit)Mr) + 5Mr,t). (229) 



From Eq.( |224D valid to first order in 6(j) we realize that rj = {(f)Q\6(f)) is such that //(t) -|- 
?7*(t) = , so that ri[t) is purely imaginary and can be reinterpreted as a change of phase 
of 00 : 

0(f, t) ~ e''W0o(r ) + 50±(f, t). (230) 

One then sees that this change of phase has no consequence on the one-body density 
matrix of the condensate, up to first order in 6(f) : 

^ |0O)(0O| + |0O)(50±| + |50±)(0O|. (231) 
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After a little algebra we turn Eg. (P26|) into a closed equation for 



ihdt 



{f,t) 



501 (r,t) 



+ 



s±if,t) 

-Sl{f,t) 



501 (f,t) ^ 

Introducing the projection operators orthogonally to 0o and (p^ 



Q 

Q* 



l-|0o)(0o| 



we have S± = QS and 



HGP + gNoQ\<Po\^Q 



(232) 



(233) 
(234) 



where the Gross-Pitaevskii Hamiltonian Hop is defined in Eg. ( P28| ). 
operator C is diagonalizable. 



(235) 

In general the 



6.1.3 Spectral properties of C and dynamical stability 

Consider an eigenvector of C with eigenvalue : 




(236) 



The free evolution of this mode, that is for 6U = , is given by the phase factor 
exp{—iekt/h) . This factor remains bounded in time provided that the imaginary part 
of efc is negative, which leads to the dynamical stability condition 

Im(efe) < for all k. (237) 

One has the following three interesting spectral properties. 

1. el is also an eigenvalue of £ . 

2. [ j is also an eigenvector of C , with eigenvalue —el . 




is an eigenvector of with eigenvalue . 



Linearization of Gross-Pitaevskii equation 



70 



The last two of these three properties can be checked by direct substitution. They can 
be viewed more elegantly as a consequence of the symmetry properties: 

1 

= -£* (238) 

-1 

= (239) 

As we shall see this last property involving is useful to write C in diagonal form. 

The first of these properties is easy to prove when 0o is real. In this case the operator 
£ is a real operator (as Q , Uq , (pl , A are real) so that complex eigenvalues come by 
pairs of complex conjugates. When 0o is complex one can use the following mathematical 
fact: if e/c is an eigenvalue of an arbitrary operator C , el is an eigenvalue of the 
operator C"^ . Here C"^ differs from £ by a unitary transform ( p39| ) so that it has 
the same spectrum as £ . This property of the spectrum is actually known in classical 
mechanics for a linearized Hamiltonian system, and the Gross-Pitaevskii equation can be 
viewed as a classical Hamiltonian equation for a continuous set of conjugate coordinates 
q = Re(0),p = lm(0) . 

As a consequence of this first spectral property of C the dynamical stability condition 
Eq.( ^37| ) can be reformulated as 

Im(efc) = for all k (240) 

that is all the eigenvalues of C have to be real to have a dynamically stable condensate 
wavefunction. We assume that this property is satisfied in the remaining part of this 
subsection |6TI| . 

6.1.4 Diagonalization of C 

As C is not a Hermitian operator the eigenbasis of C is not orthogonal (see the minus 
sign in the second line of C , due to Bose statistics; one does not have this sign in the 
BCS theory for interacting fermions). 

To write C in diagonal form the knowledge of the eigenvectors is then a priori not 
sufficient. One generally proceeds as follows: 
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Reminder: Let M be a diagonalizable but not necessarily Hermitian operator. Then 
the diagonal form of M can be written as 

M = J2mk\^k){i^k\ (241) 

k 

where is a right eigenvector of M with eigenvalue rrik : 

M\^P^) = m,\^p^) (242) 
and {ip^l is a left eigenvector of M with eigenvalue nik : 

{^PilM = m,{^Pi\ (243) 

or equivalently 

M^\^^)=ml\ij^). (244) 
The normalization of the left and right eigenvectors is such that 

= (245) 

\ipj^) is then called the adjoint vector of \ipf) . 

We apply this reminder to C . We liave already defined the right eigenvector: 

l^f ) = ( '"'=1 I , eigenvalue of £ : efc. (246) 

From the third of the above spectral properties of C we can easily obtain the correspond- 
ing left eigenvector up to a normalization factor: 

l^fc ) = I ^'',^\ ) , eigenvalue of : e* = e^. (247) 

The normalization condition Eg . (|245|) imposes 

= 1 = m{uk\uk) - {Vk\vk)]. (248) 

It is then natural to normalize the right eigenvectors in such a way that the quantity 
between square brackets is ±1 , leading to Mk = ±1 • 

We therefore group the eigenvectors of C in three families: 
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• the + family, such that {ukluf^) — {vk\vk) = +1 

• the — family, such that {uk\uk) — {vklvk) = —1 

• the family, such that {uk\uk) — {vk\vk) = . 

From the spectral property number 2 we see that there is a duality between the + 
family and the — family. Conventionally we will refer to eigenvectors of the + family 
as {uk,Vk) (eigenvalue ) and the eigenvectors of the — family will be expressed as 
{vl,u*k) (eigenvalue -e^ ). 

Generally there are only the following two members in the family: 



(:°) 


and 


( <Po ) 









(249) 



One can check that these two vectors are eigenvectors of C with the eigenvalue zero. [In 
the case of Cgp one finds in general only one zero-energy eigenmode, the missing one 
leads to the non-diagonalizability]. In general these two vectors span the whole family. 
As they are also eigenmodes of the operator with the eigenvalue zero they are actually 
their own conjugate vectors! From Eq. (|245|) we then get the important property: 

{4>o\tJ-k) = (</'ol^fc) = for all k in + family. (250) 

It is important to note that the + in the denomination " + family" refers a priori to 
the sign of {uk\uk) — {vk\vk) and not to the sign of ! 



6.1.5 General solution of the linearized problem 



We expand the unknown column vector of Eq.( ^32| ) onto eigenmodes of C . We assume 
that the only modes of the family are the ones of Eq. (|249|) ; these modes do not 
contribute to the expansion as {(t)o\S(j)±) = {4>q\S(I)]_) = . We then get the expansion 

with the coefficients 

bkit) = {{Uk\, -{vkD ( j^^^l^ll ) = / rf¥ [ul{f)SMr,t) - vl{f)S<Pl{f,t)] . (252) 
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As the second component of the expanded column vector is the complex conjugate of the 
first component the amplitudes on the — family modes are the complex conjugates of 
the amplitudes bk on the + family modes, that is bl . 

Similarly one expands the source term of Eg. (|232|) on the eigenmodes of C . The 
components on the + family modes are given by 

Skit) = J rfV [ul{f)S^{f,t)+vl{f)Sl{r,t)] . (253) 

Note the absence of — sign here, due to the fact that the second component of the source 
column vector is the opposite of the complex conjugate of the first component. Finally 
the projection of Eq.( |232| ) on the eigenmodes of the + family gives the set of equations: 

thj^bkit) = eMt) + Sk{t) (254) 
which can be integrated including the initial condition (5</)±(to) = : 

bk{t) = / ^Sk{t - r)e--'=-/^ (255) 
Jo in 

6.1.6 Link between eigenmodes of Cqp and eigenmodes of C 

The linear operators C and Cqp describe the same physical problem, so that one expects 
in particular that their spectra, which correspond to the linear response frequencies of the 
condensate, are the same. 

This expectation is confirmed by the simple result, that one can check by direct 
substitution: if {\u'^^), \v'^^)) is an eigenvector of Cqp with the eigenvalue then 
{Q\u^^) ,Q*\v'^^)) is an eigenvector of C with the same eigenvalue eu ■ 



6.2 Examples of dynamical instabilities 

We consider simple stationary solutions of the Gross-Pitaevskii equations that are dynam- 
ically unstable, that is with non-real eigenfrequencies ek/h . The situations considered 
correspond to a gas trapped in a cubic box with periodic boundary conditions; analyti- 
cal calculations can then be performed. The conclusions remain qualitatively correct for 
harmonic traps. 
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6.2.1 Condensate in a box 



The atoms arc trapped in a cubic box of size L , and we assume periodic boundary 
conditions. An obvious solution of the time independent Gross-Pitaevskii equation is 
then the plane wave with vanishing momentum, 



0o(r) 



It has a chemical potential 



II = gNo\(j)o\'^ = pog 



(256) 
(257) 



where Pq = Nq/L^ is the density of condensate atoms. 

To obtain the hnear response frequencies of the condensates we calculate the spectrum 
of Cgp ) this operator takes here the very simple form: 



J^GP 



\ 



—A + Pog 
2m 

-Pog 



Pog 



2m 



^ + Pog 



\ 



I 



(258) 



Using the translational invariance of this operator we seek its eigenvectors in the form of 
plane waves: 



(259) 











L3/2 





Within the subspace of plane waves with wave vector k , Cqp is represented by the 2x2 
non-Hermitian matrix: 



2 1,2 



CGp[k] = 



2m 



+ Pog 



\ 



-pog 



Pog 

rt'k^ 



2m 



+ Pog 



\ 



J 



(260) 



For /c 7^ this matrix can be diagonalized, giving one eigenvector of + family, with the 
eigenvalue 



2m \ 2m 



+ '^Pog 



(261) 
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and one eigenvector of 
family can be chosen as 



family with the eigenvalue —er . The eigenvector of the + 



with the correct normalization f/3 

k 




-1/4 



(262) 



V 2m 



The spectrum Eg. ( P61| ) is the so-called Bogoliubov spectrum, as it was first derived 
by Bogoliubov. Physically it is a very important result. In the limiting case of the ideal 
Rose gas { g = 0) the spectrum is the usual parabola; for g the spectrum is very 
different, and this deserves a more detailed discussion 



• Bogoliubov spectrum for g > 

In the case of effective repulsive interactions the Bogoliubov spectrum strongly differs 
from the one of a free particle in the low momenta domain Ti^k"^ /2m <^ 2pofi' as it scales 
linearly with k : 

(263) 



hk 



pog 



m 



This linear dispersion relation leads to a propagation of low energy excitations in the 
condensate in the form of sound waves with a sound velocity Cg given by 



duJi 



1 de. 



dk h dk 

or equivalently by the relativistic type formula 



Po9 
m 



(264) 



(265) 



Superfluidity is an important consequence of this linear behavior of the spectrum at 
low k , as shown by an argument due to Landau and that we explain briefly Consider a 
particle (of mass M ) sent in the atomic gas with an initial velocity u . The motion of this 
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particle can be damped by interaction with the condensate only if the particle can create 
some excitation of the condensate. Imagine that such an excitation is produced, with 
momentum k ; the particle experiences a momentum recoil of —hk and conservation of 
energy imposes 



Mu-hk 



2"'" 2M 
The velocity u has then to satisfy the inequality: 



27.2 



= Tik-u — . 266 

2M ^ ' 



u > 



k ■ u 



> ^ > C. (267) 
k 



So a particle with an incoming velocity smaller than the sound velocity can move through 
the condensate without damping, only scattering on thermal excitations of the gas can 



damp its motion. This prediction has received an experimental confirmation at MIT |^ 

At high momenta ( Tt^k"^ /2m 3> pog ) corresponding to a velocity hk/m much larger 
than the sound velocity the Bogoliubov spectrum reduces to a shifted parabola 

2m 2m 



- -iTZr + Po9 = -iTZr + ^Po9 - P- (268) 



This approximate form can be obtained by a series expansion of the general formula 



Eq.( p61| ). It can also be derived more instructively from the observation that the off- 
diagonal coupling pqq between the f/g component and the component in the 2x2 
matrix Eq. (|260|) becomes very non- resonant at high k (because the diagonal terms for 
and Vj: have opposite signs); neglecting this coupling one recovers Eq.( |268| ) with 
f/,^ ^ 1 , ^ . 

This last high energy property applies also in a non-uniform trapping potential: ne- 
glecting the off-diagonal coupling between Uk and Vk one approximates the high energy 
part of the Bogoliubov spectrum by the eigenvalues of 

-— A + [/o + 2(7iVo|0o|'-/i (269) 
2m 

which (up to the shift — /i ) is the Hartree-Fock Hamiltonian for non-condensed particles 
Eq.( |116| ) in a regime of an almost pure condensate where the density of non-condensed 
particles is negligible as compared to the condensate density Ai'ol^oP • 

From this we expect that the Hartree Fock approach is invalid for the low energy 
fraction of the non-condensed gas (energy typically less than p ); this may become a 
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problem at temperatures fc^T < /i , where one has to use the more precise Bogohubov 
approach of §0. 

• Case of a negative g 

In the case of effective attractive interactions between particles the dynamical stability 
condition Ini(e^) = is satisfied if and only if 

+ 2gpo > for all A; > 0. (270) 

2m 

If one considers the thermodynamical limit of an infinite number of condensate atoms 
with a fixed mean density po = Nq/L^ the stability condition cannot be satisfied as k 
can be arbitrarily close to in an infinite box. One may then be tempted to conclude 
that condensates with effective attractive interactions cannot be obtained experimentally, 
attractive interactions leading to a spatial collapse of the gas. 

Experiments with atomic gases can deal however with small number of atoms and 
the simplifying assumption of a thermodynamical limit is not necessarily a good approx- 
imation. In the cubix box of size L with periodic boundary conditions the compo- 
nents of the wavevector k of an atom are integer multiples of 2tt/L . The smallest but 
non zero modulus of wavevector that can be achieved is therefore 27t/L (by taking e.g. 
kx = 27T/L,ky = kz = 0). Dynamical stability condition Eq.( p7CI| ) can then be rewritten 
as 

or equivalently in terms of the scattering length as 

^ < -. (272) 

Condensates for a < can contain a limited number of atoms proportional to the size 
of the condensate. 

Condition Eq. (|271|) has a clear physical interpretation: the energy gap in the spectrum 



of a particle in the box between the ground state and the first excited states should 
be larger than the mean field energy per particle: stabilization against collapse is thus 
provided by the discrete spectrum of the atoms in the trapping potential. This condition 
can thus be qualitatively extended to the case of an isotropic harmonic trap, hu > 
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|(7|A^o/'^ho where u is the oscillation frequency of the atoms in the trap and = 
{h/muoy/'^ is the typical spatial extension of the ground state of the trap. One then 



recovers up to a numerical factor the results of § p.2.1 
6.2.2 Demixing instability 

We consider here atoms with two internal states a, h ; this model is relevant for exper- 
iments performed at JILA on binary mixtures of Rb condensates, and also (if one 
includes a third atomic internal level) experiments at MIT in Ketterle's group on spinor 
Na condensates. 

To describe the elastic interactions between the atoms with two internal states one 
needs three coupling constants, all positive in the case of ®^ Rb: Qaa and giji, for inter- 
actions between atoms in the same internal state, gab for interactions between atoms in 
different internal states: 



9aa 


a -\- a - 


a + a 


9bb ■ 


b + b- 


^b + b 


9ab '■ 


a + b - 


-> a + 6. 



(273) 



In the JILA experiment internal states a and b correspond to different hyperfine levels of 
the atoms so that inelastic collisions such as a + a a + b are either strongly endothermic 
(and do not take place) or strongly exothermic (and result in losses of atoms from the 
trap); we neglect these inelastic processes. 

Omitting for simplicity the regularizing operator in the pseudo-potential we write the 
interaction Hamiltonian between the atoms in second quantized form as 

Oaa ?t/^\ ?t/^\ 7 7 I 9bb 



(Tf 



YV'I(^)V^l(^)V^a(r)V^a(r 
+ 9ab'^Pl{r)i)l{r)i)a{f)i)i,{r\ 



(274) 



where ipa and ip}, are the atomic field operators for atoms in state a and b respectively. 
Note the absence of factor 1/2 in the a — b interaction term, which is best understood 
in first quantization point of view: all the pairs of atoms of the form a : i,b : j , where 
i running from 1 to Na labels the atoms in a and j running from 1 to Nb labels the 
atoms in b , are different so that there is no double counting of these interaction terms 
and no factor 1/2 is required. 
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Using the trick of § |5.1.3| we can rapidly derive the Gross-Pitaevskii equations for the 
condensate wavef unctions (pa in state a and (pb in state b , both wavefunctions being 
normahzed to unity. One simply has to write the Heisenberg equations of motion for the 
field operators and perform the substitution 



^6 



(275) 
(276) 



where Na^ are the number of particles in condensates a,b . As in § |6.2.1| we restrict to 
the case of atoms trapped in a cubix box of size L with periodic boundary conditions. 
We obtain the coupled time dependent Gross-Pitaevskii equations 



ihdi 



Wa 



ihdr 



-—A + Nagaa\^a\^ + A^fe^abl^fe 

2m 

-—A + Nkgbb\<^)b? + Nagab\(l>a\ 
2m 



■>b- 



(277) 



Consider now a steady state solution of these coupled Gross-Pitaevskii equations. As 
we have not introduced any coherent coupling between the internal states a and b (no 
ipb'i'a term in the Hamiltonian) (pa and (pb can have in steady state time dependent 
phase factors evolving with different frequencies: 



0a(r,t) 
(pb{r,t) 



iiafiyr )e 



(278) 
(279) 



From a more thermodynamical perspective we can also observe that the number of parti- 
cles Na and Nb are separately conserved by the three elastic interactions of Eq.( |273| ) so 
that two distinct chemical potentials fia and fib are required to describe the equilibrium 
state. 

The time independent Gross-Pitaevskii equations for and (pb.o in the box have 
the natural solutions 

(f>a,oir) = (Pb,oin = (280) 
leading to the following expressions for the chemical potentials: 

yWa = Pa,09aa + PbfiQbb (281) 
= PbfiQbb + PafiQab (282) 
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where pa^b are the condensate densities in a,b . Although we assume here that all the 
coupling constants are positive it is physically intuitive that these spatially uniform solu- 
tions should become instable when the interactions between a and b are very repulsive; 
one then feels that the two condensates a and b have a tendency to spatially separate. 

To test this expectation we linearize the time dependent coupled Gross-Pitaevskii 
equations around the steady state, setting 



Mr,t) = e-^'^''*/"[0,,o(r-') + 50,(f,t)]. 



We obtain 



ihdfi 



2m 



a,09a 



K] + Pb,Ogab[S(l)b + 



(283) 
(284) 



(285) 



and a similar equation for 6(j)b exchanging the role of a and b indices. We look for eigen- 
modes of these linear equations, with eigenfrequency e/h and a well defined wavevector 
k . This amounts to performing the substitutions 



Vae 



i{k-f—et/h) 
i{k-r—et/h) 



and equivalent changes for the b components. This leads to the eigensystem 



e — 



2m 



2m 



Va = Pa,09aa [Ua + Va] + PhfiQab [Ub + Vb] 



(286) 



and similar equations obtained by exchanging the indices a and b . Taking as new 
variables the sums and the differences between u and v and using the first identity in 
Eg. ( p86| ) we eliminate the differences as functions of the sums: 



Ua - Va 



2me 



We get from the second equality in Eq. (p86D (and a 

sums Ua + Va , Ub + Vb 



[Ua + Va] . (287) 

> b) a. two by two system for the 



2m 



/ 2me \ ' 



Id + M 



Ua + Va 
Ub + Vb 







(288) 
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where we have introduced the two by two matrix 



M 



Pafidaa PhflQab 
Pa,ogab Pbflgbb 



(289) 



This leads to the following condition for the spectrum: 



2m \ 2m 



+ 2r] 



1,2 



(290) 



where 771,2 are the eigenvalues of M . 

In the thermodynamical limit the mixture of condensates with uniform densities is 
dynamically stable provided that both eigenvalues 771,2 are positive. This is equivalent 
to the requirement that the symmetric matrix M is positive. As Qaa > here this is 
ensured provided that the determinant of M is positive: 



detM 



PafiPbfi 



QaaQbb 



2 
9ab 



> 0. 



The mixture of spatially uniform condensates is therefore stable if 

9ab < {gaagbbf^- 



(291) 



(292) 



In this case one can check that the spectrum of the + family is given by the positive 
solutions e to Eq. ( P^ ) . The spectrum of the binary mixtures of condensates is then made 
of two branches, which are both linear at low momenta with sound velocities (771,2/""^)^^^ • 
What happens when this stability condition is not satisfied ? The condensates a 



and h have a tendency to separate spatially. This happens in the JILA experiment |^ . 
Actually our model in a box is too crude to be applied to the experimental case of particles 
in a harmonic trap, the stability condition Eq.( |292| ) being marginally satisfied for ^'^ Rb; 
a numerical solution of the time dependent Gross-Pitaevskii equations is required in this 



case 



43 



The occurrence of demixing when Eq. (|292|) is violated can be connected with the 
following simple energy argument. Consider a demixed configuration with all the Na 
atoms in the left part of the box in a volume vLi^ and all the Ni, atoms in the right part 
of the box in the complementary volume {1 — v)L^ . The condensate densities vanish on 
a scale on the order of the healing length ^ of the gas, this leads to "surface" kinetic and 
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interaction energies negligible in the thermodynamical limit as compared to the volume 
interaction energy 

Nigbb 



+ 



2z/L3 2(1 - v)L^' 
We minimize this energy over v to obtain 

1 



demix 



2L3 



1/2 



(293) 



(294) 



We find that the demixed configuration has an energy lower than the one of the spatially 
uniform configuration 

1 



E, 



unif 



2L3 



K9aa + N^Qbh + 2Ar,iVfc(7, 



ab 



(295) 



precisely when the stability condition Eq. (|292|) of the uniform configuration is violated. 



6.3 Linear response in the classical hydrodynamic approxima- 
tion 

We consider in this subsection the case of a condensate in a harmonic trap. The eigen- 
modes of the linearized Gross-Pitaevskii equation can then in general be determined only 
numerically. In the Thomas-Fermi regime however approximate results can be obtained 
for the low energy eigenmodes of the system from the classical hydrodynamic approach, 
as we shall see now. 



6.3.1 Linearized classical hydrodynamic equations 

The classical hydrodynamic equations for the position dependent condensate density p(r ) 
and velocity field v{r) have been derived in §5.3.4| : 

dtp + diY[pv] = (296) 
m (j)t + V ■ giad^ V = —gTad[U + pg]. (297) 

We hnearize these equations around their steady state solution with density po and 
vanishing velocity field vq = in the unperturbed trapping potential Uq . Writing the 
trap potential perturbation of Uq : 



U{r,t) = Uo{f) + SU{f,t) 



(298) 
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and splitting p and v as 

p(r, t) 
v{r, t) 

we obtain the linearized equations: 

dt5p+ div [po5v] 
mdtSv + greid[6pg] 



po{r) + Sp{r,t) 
+ 5v{r, i) 







(299) 
(300) 



— grad^t/. 



(301) 



Taking the time derivative of the first equation we obtain a term dt5v that we can 
eliminate with the second equation. This results in a closed equation for the perturbation 
of density: 



df5p — div 



Po9 

m 



grad 5p 



div 



Po 



m 



grad 5U 



(302) 



The homogeneous part of this equation can be rewritten in a more suggestive way by 
introducing the position dependent velocity Cg given by 



mctif 



Po{r)g. 

The homogeneous part of Eg . (13021) then reads 

d^6p — div c1{f)graLd5p 



(303) 



(304) 



which corresponds to the propagation of sound waves with a position dependent sound 
velocity Cs{f) . Note that the expression ( p03| ) could be expected from the result Eq. (|265|) 
obtained in the spatially homogeneous case. 

The propagation of sound waves in a cigar shaped trapped condensate has been ob- 
served in Ketterle's group at MIT; the condensate was excited mechanically by the dipole 
force induced by a far detuned laser beam focused at the center of the trap. It is instruc- 
tive to note that to predict theoretically the velocity of sound obtained in the experiment 
one has to carefully solve Eq.( p02| ), rather than to take naively the sound velocity on the 
axis of the trap; the naive expectation would be wrong by a factor of \/2 [HI, 



6.3.2 Validity condition of the linearized classical hydrodynamic equations 

The classical hydrodynamic equations have been obtained from the time dependent Gross- 
Pitaevskii equation in hydrodynamic point of view by neglecting the quantum pressure 
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term (see §5.3.4| ). If we keep this term and linearize the resulting equation we get extra 



terms with respect to Eg. ( |301|) . One of the extra terms is 

-ASp (305) 



mpo 

which should be small as compared to the "classical" pressure term g6p : 



mpo' 



\ASp\ <t:g\6p\. (306) 



If we denote by A; a typical wavevector for the spatial variation of 6p , A6p ~ —k'^6p 
and the condition for neglecting the quantum pressure term reads 

^ « gpo- (307) 
m 

This condition can be rewritten in a variety of ways. It claims that the wavevector k 
should satisfy 

< 1 (308) 

where ^ = h/{2mpogY^'^ is the healing length of the condensate. Eg. ( |302D cannot be 
used to describe perturbations of the condensate at a length scale on the order of ^ or 
smaller. 

The validity condition can also be written as 

hk -C mCs or hkcg -C mcl = pog ~ p (309) 

In terms of the Bogoliubov spectrum for the homogeneous condensate this means that 
the wavevector k has to be in the linear part of the excitation spectrum. The energy of 
the corresponding eigenmode ~ hkcg has to be much smaller than p . Therefore only 
the eigenmodes of Eg. ( |302|) with eigenenergy much less than p are relevant: 



e < p (310) 

the higher energy modes are not an acceptable approximation of the exact eigenmodes of 
the condensate. Note that as shown in Eg. (|310|) is a necessary but not sufficient 

condition in a trap. 
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6.3.3 Approximate spectrum in a harmonic trap 

We look for eigenmodes of the homogeneous part of Eq.( p02| ) with eigenenergy e . They 
solve the eigenvalue equation 

— e^5p = div c^(f )grad(5p . (311) 

In the present Thomas-Fermi regime is a quadratic function of the coordinates as 
it is proportional to the condensate density po • We can therefore solve the eigenvalue 
equation using an ansatz for 6p polynomial in the spatial coordinates x,y,z . If 6p is 
a polynomial of total degree n , grad 6p is a polynomial of total degree n — 1 ; after 
multiplication by and action of the div operator we get a polynomial of total degree 
(n — l) + 2 — l = n. The subspace of polynomials of degree < n is therefore stable. 

For low values of the total degree n an analytical calculation is possible. For example 
in a general harmonic trap with atomic oscillation frequencies ujx,ujy,ujz one can check 
that the polynomials Sp = x,y,z (respectively) are eigenvectors with eigenvalues e = 
hujx,hu!y,hujz (respectively). These three modes correspond to the oscillation of the 
center of mass of the gas, which is exactly decoupled from all the relative coordinates of 
the particles in a harmonic trap; for this specific example the frequencies (but not the 
modes!) predicted by classical hydrodynamics are exact. These "sloshing" modes are 
used experimentally to determine accurately the trap frequencies ujx,y,z ■ 

Another important example is the case of a total degree n = 2 . One can check that 
the subclass of polynomials involving the monomials , , z"^ and 1 is stable, which 
corresponds to the ansatz 

Sp{f) = B+ Yl Ao^rl (312) 

a=x,y,z 

By inserting this ansatz into Eq. ( |311|) we arrive at the eigenvalue system 

ie/hfAo, = 2ijlA^ + u;lYAp- (313) 

This eigenvalue system can be obtained more directly as in Eq. (|216|) by a linearization 
of the equations Eq.( |212| ) for the scaling parameters Aq, around their steady state value 



Aq, = 1 . Physically the modes identified by the ansatz ( pi2| ) are therefore breathing 
modes of the condensate. The frequency of one of these modes (the one of § |5.4.3| ) has 
been measured with high precision at MIT in a cigar-shaped trap, it differs from the 
Thomas-Fermi prediction ~ (5/2Y^'^uJz (see Eq.( |218| )) by less than one percent |48[ . 
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In the case of an isotropic harmonic trap all the eigenenergies e can be calculated 
analytically (keeping in mind that modes with e > fi are not properly described by 
classical hydrodynamics !). As shown in one uses the rotational symmetry of the 
problem, as in the case of Schrodinger's equation for the hydrogen atom, with the ansatz 

6p{f) = Yr{e,<pyPiAr) (314) 

where 6, are the polar and azimuthal angles of spherical coordinates, Y/^ is the spher- 
ical harmonic of angular momentum / . The last factor Pi,n{r) is a polynomial of degree 
n in r . As P;,„ is a polynomial the recurrence relation obtained from Eq.( |311| ) for 
the coefficients of the monomials should terminate at j = n . This leads to the 
eigenfrequencies Q = e/h such that 

= (2r? + 2nl + + for any / > (315) 

where lo is the oscillation frequency of the atoms in the trap. For / = l,n = we 
recover the sloshing modes with frequency = lo . 



7 Bogoliubov approach and thermodynamical stabil- 
ity 

Imagine that we have already checked the dynamical stability of a steady state solution 
00 of the Gross-Pitaevskii equation for the condensate wavefunction. We cannot relax 
yet and be sure that the condensate will remain in state 0o in the long run. 

What is missing is a check that interaction of the condensate with the non-condensed 
cloud does not induce an evolution of the condensate far from the predicted 0o • We 
have to check what is called thermodynamical stability of the condensate as it involves 
the "thermal" , non-condensed component of the gas. 

This check will be performed in the low temperature domain ( T -C Tc ) using the 
Bogoliubov approach. We will then present examples of thermodynamical instabilities. 

We give here a summarized account of the U{1) -symmetry preserving Bogoliubov ap- 



proach developed in |^0|, In contrast to the almost general attitude in the literature 
this approach does not assume a symmetry breaking state with (■?/') ^ but consid- 
ers instead a state with a fixed total number of particles. A different U{1) -symmetry 



preserving Bogoliubov approach was developed long ago in ||52|| . 
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This allows to eliminate a technical problem of the symmetry breaking approach: if 
{ip){t = 0) 7^ the state of the system necessarily involves a coherent superposition of 
states with different total number of particles; such a state cannot be stationary (as states 
with different number of particles have also different energies) and it experiences a phase 
collapse {'ip){t) —>■ making the description of the evolution of the system more involved. 

7.1 Small parameter of the theory 

We restrict in this section to a steady state regime where most of the atoms of the gas 
are in the condensate. We split the atomic field operator as 

^(^) = 4>o{r)d^o + 5^(f) (316) 

where 0o is the condensate wavefunction and annihilates a particle in the mode 0o ■ 
The idea is to treat 5ip as a perturbation with respect to ^oo^o : Let us compare indeed 
the typical matrix elements of these two operators: 

- {SiP^Si^y/^ ~ (p')'^' (317) 
where p' is the density of non-condensed particles whereas 

00000 ~ -^o'^Vo ~ Po^ (318) 
where po is the condensate density. We will therefore assume 

p' < Po (319) 

and even more 

N' ^ J £fp'{f) <^No^N (320) 

where N is the total number of particles in the gas. Using these two assumptions a 
systematic expansion of the field equations in powers of the small parameter 

£ = {N'/Nof^ < 1 (321) 

can be performed. We give here a somewhat simplified presentation. 
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The following identities are important properties of Sip . First, dtp is the part of the 
atomic field operator transverse to the condensate wavefunction: 

d^f(j)o*{f)6ip{f) = 0. (322) 

As a consequence the bosonic commutation relation obeyed by Sip involves matrix ele- 
ments of the projector Q orthogonal to (po rather than the identity operator: 

[6^{n ), #^(r2 )] = (ri \Q\f2) = 6{n - r 2 ) - Mn Wo{r2). (323) 

Second, there should be no coherence in the one-body density matrix between the con- 
densate and the non-condensed modes, or equivalently 0o should be an eigenstate of the 
one-body density matrix (with eigenvalue Nq ): 

{aiS^P) = 0. (324) 

This last identity is used to calculate the condensate wavefunction order by order in the 



small parameter [BG, 51 . 



7.2 Zeroth order in e : Gross-Pitaevskii equation 

To zeroth order in e all the particles are supposed to be in the condensate. In steady 
state one then recovers the time independent Gross-Pitaevskii equation with the further 
approximation Nq c:^ N : 



-^A + U + gN\<Po\' 
2m 



(325) 



7.3 Next order in e : linear dynamics of non-condensed particles 

Calculation to first order in e corresponds to a linearization of the Heisenberg field equa- 
tions around 00^^00 keeping terms up to first order in 6ip . Equivalently it corresponds 
to a quadratization of the Hamiltonian around (poa^pQ keeping terms up to second order 
in Sip . 

We use this quadratization approach here. At this order of the calculation the reg- 
ularizing operator of the pseudo-potential can be neglected, divergences due the use of 
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the non-regularized g6 potential coming at the next order . We therefore take the 
Hamiltonian 



H = (fr 



(326) 



The one-body Hamiltonian hi contains the kinetic energy and the trapping potential 
energy: 

h = -—A + U. (327) 
2m 

It does not contain any — /i term as we use here in the canonical rather than grand 
canonical point of view, the total number of particles being fixed to . 

We substitute expansion (|316|) for Eq.( p26| ) and we keep terms up to quadratic in 6iIj . 
• The contribution of hi is quadratic in ■?/' so that all the terms should be kept. One 
of the contributions is 

{J £f(j)o*hi(j)o)al^a^g. (328) 

One can use the following trick to express this quantity in terms of Sip : from Eq.( |322| ) 
one sees that the total number of particles operator can be written as 



N = J d^fip^ip = + 6N (329) 
that is the sum of the number operator of condensate particles: 

= 4o^<^o (330) 
and the number operator of non-condensed particles: 

6N = [ d^f6ip^6i). (331) 



We therefore obtain 

n^^=N - SN. (332) 
• The expansion of the interaction term gives the following up to quadratic contributions: 

+ ^[<f>o* <l>o*al^al^Sip Sip + h.c.] 

+ 2g(j)o*(f)oal^5ip^Sipd^^. (333) 
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The first line of tliis expression is transformed using Eq.( |332| ): 

4o4o«</'o«</.o = (^*o - 1) = N{N - 1) - 2N6N + ... (334) 

with ~ — 1 . When we group the above term in SN with the one coming from hi 
and we replace iV by we obtain 



-6N 



(ff(t>o*hi(j)o + Ng\(t>o\ 



-^i6N (335) 



as (/)o solves the Gross-Pitaevskii equation (|325|) . This is amusing: we obtain formally a 
grand canonical Hamiltonian for the non-condensed particles, the reservoir being formed 
by the condensate particles! In the second line of Eq.( |333|) we replace a^^^d^^ by as 
the corrective term 6N would lead to a cubic contribution in 6ip . 

Another important transformation is performed by collecting the terms linear in S'ip 
from Eq.( |333|) and from the contribution of hi , leading to 

J rf¥0o*4o[^i + 9\<l>o\^N]6^ + h-c. (336) 
As 00 solves the Gross-Pitaevskii equation, the operator between brackets, when acting 
on the left on 0o* , gives a contribution /i0o* orthogonal to Sip (see Eg . (13221) ). The 
sum of all the terms linear in 6'ip therefore vanishes! We shall see later that this could 
be expected from Eq. (|324|) . 



7.4 Bogoliubov Hamiltonian 

We now collect all the terms of H up to quadratic in S'ljj and express them in terms of 
the field operator 

A(r ) = ^dl6^if)- (337) 

This operator commutes with the operator N giving the total number of particles: it 
transfers one non-condensed particle into the condensate. As the operator , it is 
transverse to 0o (see Eq.( |322| )). By definition of the condensate wavefunction A has a 
zero mean (see Eq. (|324|) ). 

In general it is difficult to exactly eliminate 6tp in terms of A . Fortunately at the 
present order of the calculation we can use the assumption of a very small non-condensed 
fraction so that one has for example 

6tp^6tp ~ 6tp^a^,N~^d\6'^ = . (338) 
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To the same order of approximation A obeys the same commutation relation Eq. (|323|) 
as 5?/' . 

The final result can be written in term of the operator C introduced in §|^, with the 
approximation Nq c:^ N : 



^quad = f{N) + \j d'f{A\ ~A)C j 



(339) 



where the function / is specified in [pT[] . 



From this quadratic Hamiltonian the Heisenberg equations of motion for the field A 
have the suggestive form 

'''' it)- (3«) 

We note that an hypothetic term of i/quad linear in A would give rise to a source term 
in Eq . (|340|) preventing one from satisfying Eq . (13241) at all times! 




The result Eq. (|34CI|) is really a crucial one. It shows that the linearized evolution of the 



non-condensed part 6tp of the atomic field is formally equivalent to the linearized response 
of the condensate to a classical perturbation (e.g. of the trapping potential) derived from 
the Gross-Pitaevskii equation; both treatments indeed exhibit the same operator C (see 
Eq.dH)). 



All the machinery of can therefore be used. We expand the field operator A on 
the eigenmodes of C . We assume here dynamical stability and that the only eigenmodes 
in the family are the zero-energy modes ((/)o, 0) and (0, (j)o*) to which the field A is 
orthogonal according to Eq. (|322|) . We therefore get an expansion similar to Eq. (|251D : 




= E h \+bl ^f/ (341) 

with the important difference that the coefficients in the expansion are now operators: 



ul{r)k{f)-vl{r)k\r)\. (342) 



From the normalization condition between the eigenvectors of C and their adjoint vectors 
(see §|]) one shows that the 6^ obey bosonic commutation relations: 

[hk, h\,] = 5k,k' and [h, h'] = 0. (343) 
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bk corresponds formally to an annihilation operator; as \vk) in general does 
not simply annihilate a particle as it is a coherent superposition of A (which transfers 
one non-condensed particle to the condensate) and of A''^ (which transfers one condensate 
particle to the non-condensed fraction). One then says that bk annihilates a quasi-particie 
in mode k . 

Finally we rewrite the Hamiltonian Eg. ( |339|) in terms of the bk 's: 

^quad = ^o(iV) + ^ ekblbk. (344) 

fcG+family 

We recall that is the eigenenergy of the mode {uk, Vk) of the + family. Our quadratic 
Hamiltonian describes a gas of non-interacting quasi-particles: this is the so-called Bo- 
goliubov Hamiltonian. 

The ground state of -f^quad is obtained when no quasi-particle is present, it corresponds 
to all the modes k being in vacuum state: 

i^quad|0) = ^o(iV)|0) (345) 

where 

hk\Q) = Vfc. (346) 

Eq is therefore the Bogoliubov approximation for the ground state energy of the gas. 
To get a finite expression for Eq one has to include the regularizing operator in the 
pseudo-potential. 

The excited states of the system are obtained in the Bogoliubov approximation by 
successive actions of the fe^ 's. For this reason b\ is said to create an eiemeritar j excitation 
k in the system, to distinguish with collective excitations involving all the particles of the 
condensate (induced e.g. by a perturbation 5U of the trapping potential). We emphasize 
again that the elementary excitations of the gas have the same frequency ek/Ti as the 
collective excitations in the linear response domain ( 5U small enough). This intriguing 
property is valid only at the presently considered regime of an almost pure condensate 
(T<T,). 

7.5 Order : corrections to the Gross-Pitaevskii equation 

Expanding the Heisenberg field equations for ■?/; keeping terms up to N^/'^e'^ one can 
calculate the first correction to the prediction Eq.( |325| ) for the condensate wavefunction. 
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This correction includes (i) the fact that the number of condensate particles Nq rather 
than the total number of particles should appear in the Gross-Pitaevskii equation, 
and (ii) the mechanical back-action of the non-condensed particles on the condensate in 
the form of mean field potentials. 



The calculations are a bit involved |5T[ and require the use of the regularizing operator 
in the pseudo-potential (a fact realized in ||53l but not yet in |^). We give here only the 
result. The condensate wavef unction is given by an expansion 



(347) 



where (po^^^ , zeroth order approximation in e , is the solution of Eq. (|325|) . The correction 



is of order ; its component on 00*'°'* is purely imaginary (as both (pQ and ipo 



(0) 



are normalized to unity) and can be considered as a (not physically relevant) change of 
global phase of 0o''°^ ! the part of ^o^^"* orthogonal to 0o''°^ is given by: 

\ f QS \ 



Q*(Po- ^* J V -Q*s* J 

where Q projects orthogonally to (po and where the source term S is equal to 



(348) 



- {l + {6N))g\<f>o^'\f)W\r) 

+ 2g{A\r)A{f))(l)o^^\f) + g [d^ (s(A(f - s/2)A(r + s/2 

'rf^s |0o^°H5')p([At(.-')0o(°H^"') + A(s-')0o^°^*(5')l A(f)). 



(349) 



The first line of this expression contains the effect of the depletion of the condensate, 
the number of non-condensed particles {6N) being calculated in the Bogoliubov approx- 
imation, see discussion in § [7.7| . The other terms are mean field terms, among which one 
recognizes the Hartree-Fock contribution 2g{A^{r)A{f)) already obtained in §^ 



7.6 Thermal equilibrium of the gas of quasi-particles 

In the Bogoliubov approximation the quasi-particles behave as an ideal Bose gas; such a 
gas can reach thermal equilibrium only by contact with a thermostat. There is no such 
thermostat in the experiments on trapped Bose gas, relaxation to thermal equilibrium 
has to be provided instead by interactions between the atoms. 
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Fortunately the full Hamiltonian Eq.( p26| ) contains terms cubic and quartic in Sip : 
when expressed in terms of the b 's and 6^ 's they correspond to interactions between the 
quasi-particles which will provide thermalization. Two situations can then be considered, 
depending on the sign of . 

• "Good" case: > for all k in + family. We assume for simplicity thermal 
equilibrium in the canonical point of view (which should be equivalent to the micro- 
canonical point of view in the limit of large number of particles) with a N -body density 
matrix 

Pi,...,N = ^e-^^ ^e-'^^q-d. (350) 

We suppose therefore that the interactions between quasi-particles, essential to ensure 
thermalization, have a weak effect on the thermal equilibrium state. From Eq. (|350|) we 
finally obtain the mean number of quasi-particles in mode k : 

iblh) = (351) 

This good case corresponds to a thermodynamically stable condensate in state 0o • 

• "Bad" case: there is a mode ko in the -|- family such that < . In this case the 
quadratic Hamiltonian -^quad contains a harmonic oscillator of frequency ujko = Icfcol/^ 
having formally a negative mass M : 

- hofblK = ^[Pl + (352) 

where P^q and Qk^ correspond formally to a momentum and position operator. By 
collisions with the quasi-particles of positive energy the mode k^ can loose energy which 
increases its own excitation; if the number of quanta in the mode can become comparable 
to the process of thermalization of the gas may lead the condensate to a state different 
from the predicted 0o • 

This phenomenon of thermodynamical instability should not be confused with dynam- 
ical instability of the condensate, where e^p is complex; e.g. the case of a purely imaginary 

corresponds formally to an oscillator in an expelling potential, [Pko~\^ko\'^Q\^/ (2M) . 

7.7 Condensate depletion and the small parameter (pa^)^/^ 

We assume the situation of thermodynamical stability. We calculate the mean number 
of particles out of the condensate, that is in all the modes orthogonal to the con- 
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densate wavefunction 0o • 

{6N) = I £r{5%l)\r)5%^{f)). (353) 



In this way we can calculate explicitly the small parameter of the present theory given in 
Eq.dH). 

To lowest order in the Bogoliubov approximation we can replace 5il) by A in the 
above expression: 

{5N) - [ Sr{k\r)k{f)). (354) 



We replace A by its modal expansion; using the approximation in Eq.( ^5(J| ) the only 
terms with non-zero mean are ^\hk) given by Eq.( |351| ) and 

{khl) = (blh) + 1. (355) 

We therefore obtain: 

{6N) ~ E {blbk)[{uk\uk) + {vk\vk)] 

fee+family 

+ E i^k\vk). (356) 

fcG+family 

The contribution of the occupation numbers {b\bk) corresponds to the so-called thermal 
depletion of the condensate, as it is non-zero only at finite temperature. The contribution 
of the +1 coming from the identity ( |355| ) is the so-called quantum depletion of the 
condensate: it expresses the fact that, even at zero temperature, there is a finite number 
of particles out of the condensate due to atomic interactions, contrarily to the ideal Bose 
where all the Vk 's vanish. 
It is instructive to calculate the depletion of the condensate in the homogeneous case 
in the thermodynamical limit, replacing the sum over the wavevector k characterizing 
each mode of the + family by an integral. From Eq.( |262| ) we calculate 



{uk\uk) + {vk\vk) = 1 + 2{vk\vk) = (357) 



At zero temperature we obtained for the non-condensed fraction of particles the following 
integral: 

((5iV) 16 / .^1/2 



N 



9^ + 1/2 ^ 
,g(l + g2)V2 ^ 



(358) 
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where we have integrated over the sohd angle in spherical coordinates and we have made 
the change of variable 

^ = 2pgq'. (359) 
2m 

This integral can be calculated exactly: 

In this way we obtain a very important small parameter of the theory, (pa^)^^^ , which 
characterizes the domain of a weakly interacting Bose gas. This small parameter is similar 
to the condition obtained in Eq. (^) of §^ already with a totally different point of view, 
a condition of Born approximation on the pseudo-potential! In the typical experimental 
conditions of MIT for sodium atoms we have p ~ 10^^ cm ~^ and a ~ 25 A, which leads 
to a small parameter (pa^)^^"^ ~ 10~^ . The result ( |360| ) also gives us the opportunity to 
recall that the number of non-condensed particles should not be confused with the number 
of quasi-particles for an interacting Bose gas: one notes here that the second quantity 
vanishes at T = whereas the first one does not. 

At finite temperature there is, in addition to the quantum depletion, a thermal deple- 
tion of the condensate: 

(SN) (SN).^. 32 f+°° , g(g2 + i/2) 



N 



The integral over q depends only on the parameter pg/lksT) . At low temperature 
ksT <^ pg this parameter is large so that the modes with a g ~ 1 have a very low 
occupation number: one can neglect q as compared to one and one-half (which amounts 
to restricting to the linear part of the Bogoliubov spectrum) and one obtains a small 
correction to the zero temperature case: 



2pg 



+ . 



(362) 



At high temperature ksT ^ pg the major fraction of the populated Bogoliubov modes 
have a q much larger than unity, so that we can now neglect one and one-half as compared 
to (which amounts to restricting to the quadratic part of the Bogoliubov spectrum). 
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To this approximation we recover the ideal Bose gas result 

where AdB is the thermal de Broglie wavelength and ( stands for the Riemann Zeta 
function (see ^). For a given temperature indeed the density of non-condensed particles 
in presence of a condensate saturates to its maximal value C(3/2)A^b for the ideal Bose 
gas in a box (see Eq.(|^) and Eq.(|^)). 

In the high temperature regime fc^T ^ pg our Bogoliubov approach can therefore be 
valid only if pA^g ^ C(3/2) . Both inequalities on p can be satisfied simultaneously only 
if 2(^(3/2)a/AdB ^ 1 • Amusingly this condition is similar to the condition aAk <^ 1 
obtained in §^ (see Eg . (|96|) ) for the validity condition of the Born approximation for the 
pseudo-potential. 



7.8 Fluctuations in the number of condensate particles 



The Bogoliubov theory that we have presented also allows a calculation of the fluctuations 
of A'^o in the canonical ensemble. This has the advantage of removing the effect of 
fluctuations in the total number of particles present in the grand-canonical ensemble, a 
"trivial" contribution to the fluctuations of A^^o . 

As the total number of particles is fixed to A^ the variance of A^o is exactly equal 
to the variance of 6N , number of non-condensed particles. The mean value of 6N has 
been given in the previous section, we now have to calculate the mean value of its square: 



dfi J dr2 {S^'H'^i )'^^(^i )S'^\^2 )6-iJj{f2 )) 
{6N) + fdfi f df2 {6tij\fi )6i:\r2 )Sij{fi )(5V'(f2 )) 



(364) 
(365) 



where we have used the commutation relation Eq. (|323|) and the fact Eq.( p22| ) that (j) is 
orthogonal to Sip . To lowest order in the Bogoliubov approximation we can replace 6i/j 
by A in the above expression and take the approximation Eq. ( P5(]| ) for the equilibrium 
density operator. As -ffquad is quadratic in the field A , Wick's theorem can be used. We 
finally find for the variance of the number of non-condensed particles: 



Var(5iV) ~ {SN) + dn dr 



(At(fi)A(f2)) + (A(fi)A(f2)) 



(366) 
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The variance in Eg. ( |366| ) is simple to calculate in the homogeneous case of a gas in 
a cubic box of size L with periodic boundary conditions. The eigenmodes {u, v) are 
simply plane waves Eq.( p59|) with real coefficients Uk, Vk given in Eq. (|262|) and depending 
only on the modulus k of the wavevector k . We find for the two correlation functions 
appearing in Eq.( p66| ) the simple expression: 

{^\n)A{f,)) = -^E[(^' + ^.'K + ^.1e^'-^"""^ (367) 

(A(fi)A(f2)) = l^J2UkVk{l + 2nk)e'^'<''^-'''^ (368) 

where we have introduced the mean occupation numbers Uk = {blbk) ■ After spatial 
integration of the square modulus of these quantities we obtain 

Var(5iV) = {6N) + [{u^ + V,') Uk + if]' + UlV^{l + 2nkf (369) 

MO 

where the mean number of non-condensed particles {SN) is already given in Eq. (|356|) : 

m = E [Ul + Vf ) Uk + Vl (370) 



We first apply formula (|369|) to the limiting case of zero temperature. All the occu- 
pation numbers Uk vanish. For the ideal Bose gas ( (7 = ) all the Vk 's are zero and 
the variance of 5N vanishes as expected, since all the particles are in the ground state 
of the box. For the interacting Bose gas we find 

Va.(.iV)(T = 0) = ii:^I(lV^ (371) 

where q is given as function of k by Eq.( |359D . In the thermodynamical limit we replace 
the discrete sum by an integral and this leads to a variance scaling as the number of 
particles for a fixed density: 

Var(5iV)(T = 0) ^ / 3\V2 



27^^'\pay . (372) 



In the regime of validity of the Bogoliubov approach one has pa^ <^ 1 so that the 
fiuctuations of Nq are sub-poissonian. 
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The situation can be totally different at finite temperature. Consider first the case of 
the ideal Bose gas fS^ : 



Var(5iV) = Y^rikil + Uk) 

= g4sinh^U/2) ^''^^ 

with ek = h^k'^ / {2m) . In the thermodynamical limit one may be tempted to replace 
in the usual manner the sum over k by an integral. This leads however to an integral 
divergent in = : the integrand scales as , which is not compensated by the 

Jacobian k'^ of three-dimensional integration in spherical coordinates. In this case the 
contribution of the sum in the thermodynamical limit is dominated by the terms close to 
/c = where (3ek <^ 1 and the function sinh can be linearized. We then obtain 

Var(5iV),=o-(^)'E4- (374) 

In this expression we have introduced the kinetic energy difference between the ground 
state and the first excited state for a single particle in the box: 

and the sum ranges over all the vectors n with integer components and a non-vanishing 
norm n . By a numerical calculation we obtain 

E 1 = 16.53 . . . (376) 

A remarkable feature is that the variance of Nq scales as that is as the volume of the 
box to the power 4/3 , or equivalently as the number of particles to the power 4/3 in 
the thermodynamical limit. This is much larger than in the thermodynamical limit. 

Do the fluctuations of Nq remain large in presence of interactions ? As the spectrum 
efc is linear at small k the divergence of nl is only as l/k"^ , which has a finite integral in 
three dimensions. However the mode functions Uk, Vk are also diverging in A; = , each 
as l/k^/'^ , so that one recovers the 1/k^ dependence of the summand close to k = . 
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As in the ideal Bose gas case we replace in the thermodynamical limit the summand by 
its low k approximation: 

- ^ (378) 
V, - (379) 



and we keep in the summation the most diverging terms. We finally obtain |55, Bq 



Var(5iV),>o - M ^ ) T.^,- (380) 
Remarkably this result differs from the ideal Bose gas case Eq. ( |374D by a factor 1/2 only: 



fiuctuations of A^^o remain large. The fact that Eq. ( |380| ) does not depend on the strength 
g of the interaction is valid only at the thermodynamical limit and indicates that the 
limit g and the thermodynamical limit do not commute. 

7.9 A simple reformulation of thermodynamical stability condi- 
tion 

The thermodynamical stability condition simply requires that the Bogoliubov Hamilto- 
nian Hquad is the sum of a constant (function of N ) and of a positive quadratic operator 
in the field variables. In the diagonal form ( |344| ) positivity is clearly equivalent to the 
requirement positive for all k in the + family. How to express this condition from 
the non-diagonal form ( p39| )? We simply have to rewrite Eq. (|339|) as 



i^quad = fiN) + \j d'r{k\ A)vC j (381) 

where rj is the operator 



?) (382) 



so that r]C is an Hermitian operator. 
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The thermodynamical stability condition is therefore equivalent to the requirement 
that rjC is positive: 

riC > 0. (383) 

More precisely, as A is orthogonal to 0o , V'^ has to be strictly positive in the subspace 
orthogonal to (|0o))O) and (0, |0o*)) • 

We can give a simple physical interpretation of this condition: 0o has to be a local 
minimum of the Gross-Pitaevskii energy functional 



*]=N d^r 



— |gr;d0p + U{fmr)\' + ^Ng\m\' 



(384) 



which is the expression of §|] with the approximation No c:^ N . Let us consider indeed 
the variation 6E of E from Eq = E[(j)Q,(f)Q*] up to second order in a small deviation 
50 of (p from (po . 

The terms linear in 50 are given by: 



- — grad 00* ■ grad Scj) + U (f )0o*(50 + A^5'0o*^0o'^0 + cc. 
2m 



(385) 



By integration by parts and using the fact that 0o solves the Gross-Pitaevskii equation 
Eq.( p25| ) we rewrite this expression as 

5E(^) = iV/i[(0o|50) + (50|0o)]. (386) 
As both and 0o are normalized to unity 6(f) actually fulfills the identity 



(0O|50) + (50|0O) = -(50150). 

The a priori first order energy change is a posteriori of second order: 

5^(1) = ~Nfi{S(j)\S(j))\ 
The terms a priori quadratic in 50 are given by: 



(387) 



(388) 



5^(2) = Njd^ 



rad50* ■ grad 50 + f/(f )50*50 + 2A^^|0o|^50*- 



2m 



(389) 
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We transform this expression by splitting 6(f) in a part parallel to 0o and a part orthog- 
onal to 00 : 

50(f) = 700 (r) + 50±(f ). (390) 

Using integration by parts, the Gross-Pitaevskii equation and the fact that operator C 
contains projectors orthogonally to 00, 0o* we are able to write 5E^'^^ as 

5^(2) = |^|2^^ + 1 + ^*)2^2^ I ^3^-.|^^|4 ^ ^ ^*)^2^ I rfV|0o|2(0O*50± + C.C.) 

+iiV I rfV(501,50^)(r/£ + ^Id) ( ) • (391) 

Note that we had to add fi times the identity matrix Id to r]C as Eg. ( P89| ), contrarily 
to TjC , does not contain any term proportional to . From Eg .( 13871 ) we see that 7 + 7* 
is actually of second order in 6(f) ± so that it can be set to zero in Eg. ( ^91| ). 

Summing the a priori first and second order energy changes we see that 5E^^^ exactly 
cancels the terms involving explicitly n in Eg.( p91|) so that we arrive at 

5Ec^\n j d'f{6(f)l, 6(f)^)vC ( ) • (392) 

Thermodynamical stability that is positivity of rjC is therefore eguivalent to the Gross- 
Pitaevskii energy functional having a local minimum in 0o . 

7.10 Thermodynamical stability implies dynamical stability 

As we show now the positivity of rjC automatically leads to a purely real spectrum for C , 
that is to dynamical stability. Consider an eigenvector {u, v) of C with the eigenvalue 
e. Contracting the operator rjC between the ket (|u), |f )) and the bra {{u\,{v\) we get 

((n|,(i;|)r//:(^ j^l ^ = e [{u\u) - {v\v)] . (393) 

The matrix element of rjC is real positive as rjC is supposed to be a positive hermitian 
operator. We now face two possible cases for the real guantity {u\u) — {v\v) : 

• {u\u) — {v\v) = . In this case rjC has a vanishing expectation value in (|u), \v)) ; 
as rjC is positive, {\u), \v)) has to be an eigenvector of rjC with the eigenvalue 
zero; as rj is invertible we find that (|m), \v)) is an eigenvalue of C with the 
eigenvalue 0, so that e = is a real number. 
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{u\u) — {v\v) > : we get e as the ratio of two real numbers, so that e is reaL 



7.11 Examples of t her mo dynamical instability 
7.11.1 Real condensate wavefunction with a node 

Let us assume that the solution of the Gross-Pitaevskii equation is a real function 0o(^) • 
To decide if this solution is thermodynamically stable one has to check the positivity of 
the operator rjC . Consider an eigenvector of rjC with eigenvalue e : 



(394) 



This e should not be confused with the quasi-particle energies as rjC and C have 





different spectra. By performing the sum and the difference of the two lines of Eq.( 394 ) 
we get decoupled equations for the sum ipg = u + v and the difference ipd = u — v : 



^ + u{r) + Ng<Pl{r) + 2NgQ^l{f)Q - /i 
2m 



2m 



+ U{r) + Ng(t)l{r\ 



/i 



(395) 
(396) 



Both operators involved in these equations have to be positive to achieve thermodynamical 
stability. Note that for > the positivity of the second operator Eq. ( |396| ) implies the 
positivity of the first one Eq. (|395|) as gQ<t>1{f)Q is positive. 

We therefore concentrate on Eq.( ^96| ). It involves the Gross-Pitaevskii Hamiltonian 



nGP = ^ + U{r) + Ng<g{r 
2m 



(397) 



An obvious eigenvector of this Hamiltonian is ipd = 0o with eigenvalue £ = , as 0o 
solves the Gross-Pitaevskii equation! The condition of a positive e in Eq. ( |396| ) simply 
means that (pQ should be the ground state of Ti-cp ! 

We can then invoke a theorem claiming that the ground state of a potential has no 
node IQ. If 4>o{r) has a node it cannot be the ground state of Hop ■ The ground state 
of Ti-GP has therefore an eigenenergy e lower than the one of 0o , that is lower than 
zero, so that the operator t]C is not positive and there is no thermodynamical stability. 
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As an example consider in a harmonic trap with eigenaxis z , a solution of the Gross- 
Pitaevskii equation even along x and y but odd along z , so that it vanishes in the plane 
z = . Such a solution exists, for g > : within the class of real functions (p odd along 
z and even along x ,y , the Gross-Pitaevskii energy E[(f),(f)] , bounded from below, has 
a minimum, reached in 0o , and this (po then solves the Gross-Pitaevskii equation. This 
solution however is no longer a local minimum of E[(f), 0*] when one includes all possible 
deviations of (p from 0o (complex and with no well defined parity along z ). 



7.11.2 Condensate with a vortex 

Can we get a thermodynamically stable condensate wavefunction with a node? To beat 
the results of the previous subsection we now assume 0o to be complex. 

A particular class of complex wavefunctions with a node are condensate wavefunctions 
with vortices. A vortex is characterized (i) by a nodal, not necessarily straight, line in 
00 (the center of the so-called vortex core) and (ii) by the fact that the phase of 0o 
changes by 2g7r , q non-zero integer, along a closed path around the vortex core ( q is 
the so-called charge of the vortex). This second property means that the circulation of 
the local velocity field (defined in § |5.3.3| ) around the vortex core is 2Tihq/m . 

A condensate wavefunction can have several vortices; the change of the phase of 0o 
along a closed path is now 27rgsum where qsnm is the algebraic sum of the charges of the 
vortex lines enclosed by the path. 

It has been shown that a condensate wavefunction with a vortex in a harmonic trap 



is not thermodynamically stable In the limit of vanishing interaction between the 
particles ((7 = 0) this is clear indeed. Suppose that the trap is cylindrically symmet- 
ric with respect to z . |0o) can be chosen as (|n^. = l^Uy = 0, ra^ = 0) + i\n^ = 
0,ny = l,nz = 0))/-\/2 where \nx,ny,nz) is the eigenstate of the harmonic oscillator 
with quantum number along axis a { a = x,y, z ). The chemical potential is simply 
2huJx y + \%ujz where Ua is the atomic oscillation frequency along axis a . One then 
finds that {\u) = In^ = 0,ny = 0,11^ = 0), \v) = 0) is an eigenvector of rjC with the 
strictly negative energy e = —hujx,y ■ 

What happens in the opposite Thomas-Fermi regime of strong interactions? An intu- 
itive answer can be obtained in a 2D model of the Gross-Pitaevskii equation, assuming 
for simplicity a quasi-isotropic trapping potential and restricting to the following class of 



Bogoliubov approach 



105 



condensate wavef unctions: 

Mx,y) = (j)siow{x,y) td.nh[K\f - aR\]e'^^^ . (398) 

In this ansatz (f)s\ow{x, y) is the usual square root of inverted parabola Thomas-Fermi ap- 
proximation for a condensate wavefunction without vortex, with a radius R ; the tanh[ ] 
represents the correction to the modulus of 00 due to the vortex core (of adjustable 
position dR and inverse width k ) ; 6*^/2 is the polar angle of a system of Cartesian co- 
ordinates (X, Y) centered on the vortex core, and represents (approximately for a 7^ ) 
the phase of the unit-charge vortex. 

One then calculates the mean energy of 0o , with the simplification that (f)siow{x,y) 
varies very slow at the scale of , and one minimizes this energy over k . This leads to 
the inverse size of the vortex core on the order of the local healing length of the condensate: 



0.59 



/i — -muj'^{aRY 



(399) 



m 

The mean energy of (/)o (|384|) is now a function of the position of the vortex core only, 

E = ^no vortex + W (a) (400) 



where E^o vortex is the energy of the condensate with no vortex and 



(M 1 1 . 9x 



21n2 + l ^f^o . . 

h In h ln(l - a ) 

6 nuj 



(401) 



with V = 0.49312 and /io the chemical potential in the absence of vortex. This function 
W represents an effective potential seen by the vortex core. As shown in Fig.|l^ this 
potential is maximal at the center of the trap so that it is actually an expelling potential 
for the vortex core: shifting the vortex core away from the center of the trap lowers the 
condensate energy. 

A method to stabilize the vortex is to rotate the harmonic trap around 2 at a fre- 
quency Q (the trap is anisotropic in the x — y plane otherwise rotation would have 
no effect). Thermodynamical equilibrium will now be obtained in the frame rotating at 
the frequency Q , where the harmonic trap is time independent. As this frame is non 
Galilean the Hamiltonian and therefore the Gross-Pitaevskii energy functional have to 
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be supplemented by the inertial energy term —VtLz per atom, where is the angular 
momentum operator along z . The effective potential W{a) gets an extra term: 

Wn{a) = Wn=Q{d) - Nhnil - a^f (402) 

where W^=q is the result ( [40 1|) in the absence of rotation. As shown in Fig.p!^ this extra 
term can trap the vortex core at the center of the harmonic trap if VL is large enough. 

What happens if VL is increased significantly ? It becomes favorable to put more 
vorticity in the condensate. As the vortices with charge larger than one are unstable the 
way out is to create several vortices with unit charge. This can be analyzed along the 
previous lines by a generalized multi- vortex ansatz, as discussed in [^. A condensate 
with vortices has been recently obtained at the ENS in a rotating trap . 



Another philosophy was followed at JILA: rather than relying on thermal equilibrium 
in a rotating trap to produce a vortex they used a "quantum engineering " technique 
]62[| to directly induce the vortex by giving angular momentum to the atoms through 



coupling to electromagnetic fields |^. It has also been suggested to imprint the phase of 
the vortex on the condensate through a lightshift induced by a laser beam whose spatial 



intensity profile has been conveniently tailored ||6^. Such an imprinting technique has 
successfully led to the observation of dark and gray solitons in atomic condensates with 
repulsive interactions in Hannover and in the group of W. Phillips at NIST. All these 



techniques illustrate again the powerfulness of atomic physics in its ability to manipulate 
a condensate. 



8 Phase coherence properties of Bose-Einstein con- 
densates 

Consider two Bose-Einstein condensates prepared in spatially well separated traps and 
that have 'never seen each other' (e.g. one rubidium condensate at JILA and one rubidium 
condensate at ENS). It is 'natural' to assume that these two condensates do not have a 
well defined relative phase. However the trend in the literature on Bose condensates is 
to assume that the two condensates are in a coherent state with a well defined relative 
phase, the so-called 'symmetry-breaking' point of view. So imagine that one lets the two 
condensates spatially overlap. Will interference fringes appear on the resulting atomic 
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Figure 14: In a 2D model, effective potential energy W oi a vortex in a quasi-axisymmetric 
harmonic trap as function of the distance aR of the core from the trap center, for /io = 80?ia; . 
(a) = and (b) 17 = 0.045ci; . The unit of energy is Nhuj where to is the oscillation 
frequency of the atoms in the trap. 



density or not ? 

One of the goals of this chapter is to answer this question and to reconcile the symmetry 
breaking point of view with the 'natural' point of view. 



8.1 Interference between two BECs 

At MIT a double well trapping potential was obtained by superimposing a sharp barrier 
induced with laser light on top of the usual harmonic trap produced with a magnetic 
field. In this way one can produce two Bose-Einstein condensates, one on each side of the 
barrier. The height of the barrier can be made much larger than the chemical potential 
of the gas so that coupling between the two condensates via tunneling through the wall 
is very small. In this way one can consider the two condensates as independent. 

One can then switch off the barrier and magnetic trap, let the two condensates bal- 
listically expand and spatially overlap. One then measures the spatial density of the 
cloud by absorption imaging. This spatial density exhibits clearly fringes (see figure 
ISl ). These fringes have to be interference fringes, as hydrodynamic effects (such as sound 
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waves) are excluded at the very low densities of the ballistically expanded condensates. 
We show here on a simple model that we indeed expect to see interference fringes in such 
an experiment, even if the two condensates have initially no well defined relative phase. 





Absorption 



50% 



Figure 15: Interference fringes between two condensates observed at MIT [36| 



8.1.1 A very simple model 

In our simple modelization of an MIT-type interference experiment we will concentrate on 
the positions of the particles on an axis x connecting the two condensates so that we use 
a one-dimensional model enclosed in a box of size L with periodic boundary conditions. 
We assume that the system is initially in the Fock state 

\^) = \j:K,j:h) (403) 

with N/2 particles in the plane wave of momentum hka and N/2 particles in the plane 
wave of momentum hkb : 

{x\ka,b) = ^exp [ika,bx] . (404) 
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We assume that one detects the position of all the particles. What will be the outcome ? 
As the numbers of particles are exactly defined in the two modes a and b the relative 
phase between the atomic fields in the two modes is totally undefined. 

8.1.2 A trap to avoid 

If we calculate the mean density in the state given by ( [40 3|) we find a uniform result 

{'4j^{x)ip{x)) = N/L (405) 

and we may be tempted to conclude that no interference fringes will appear in the beating 
of two Fock state condensates. 

Actually this naive statement is wrong. Interference fringes appeared in a single 
realization of the experiment at MIT. We have therefore to consider the probability of the 
outcome of a particular density profile in a single realization of the measurement and not 
the average of the density profile over many realizations of the experiment. Indeed we 
will see that by interfering two independent Bose-Einstein condensates we get interference 
fringes on the density profile in each single realization of the experiment but the position 
of the interference pattern is random so that by averaging the density profile over many 
realizations we wash out the fringes. 

We wish to emphasize the following crucial point of the quantum theory: Whatever 
single-time measurement is performed on the system all the information about the out- 
comes of a single realization of the measurement procedure is contained in the N— body 
density matrix, here 

p=|^)(^|. (406) 

Indeed the only information we can get from quantum mechanics on a single realization 
outcome is its probability P , which can be obtained from p by 

P = Tr[6p] (407) 

where the operator O depends on the considered outcome. E.g. in our gedanken exper- 
iment P is the probability density of finding the N particles at positions Xi,X2, ....xn 
and the operator O is expressed in terms of the field operator as 

O = -^ij^{xi)....^P^XN)i^{xN).-i^{xi). (408) 
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In a first quantized picture this corresponds to the fact that the probabihty density P is 
equal to the modulus squared of the A^— body wavefunction. 

The complete calculation of the A^— body distribution function P{xi, . . . ,xn) for 
the state |\E') in Eq. ( |403|) is involved and we will see in the coming subsections how to 
circumvent the difficulty. But we can do a simple calculation of the pair distribution 
function of the atoms in state l^') : 



p{xi,X2) = {^\-4j\xi)ij^{x2)^{x2)ij{xi)\^) 

= mx2mx,m\\\ 



(409) 
(410) 



We expand the field operator on the two modes (f)a,b and on other arbitrary orthogonal 
modes not relevant here as they are not populated in |\E') : 



ip{x) = d{x\ka) + b{x\kb) + . . . 



(411) 



where d and b annihilate a particle in state ka and k^ respectively. We obtain 

^Pix2)'^ixi)\^) = 



+ 



y 

N 

y 

2 



A^ 

y 

A^ 

y 



1/2 jY jY 

{X2\ka){xi\ka)\— - 2 : /c^, — : kb) 



- 1 



1/2 

{x2\kb){xi\kb)\— : k 



1 A^ A^ 
{X2\ka){xi\kb) + {X2\kb){xi\ka)\\ — - 1 : ka, — - I : kb). {412) 



The last line of this expression exhibits an interference effect between two amplitudes, 
that could not appear in the previous naive reasoning on the one-body density operator 
Eq.( [405| )! In the limit A^ ^ 1 and using the fact that the populated modes are plane 
waves the pair distribution function simplifies to 



P{Xi,X2) 



1 



COS 



{ka - kb){xi - X2) 



(413) 



Tj r ■ 2 

This function exhibits oscillations around an average value equal to the square of the 
mean density. The oscillations are due to the interference effect in Eq.( f412[ ): they favor 
detections of pairs of particles with a distance \xi — X2I equal to 2 mr / \ ka — kb\ { ti 
integer) and they rarefy detections of pairs of particles with a distance {2n+l)7i/\ka — kb\ . 
We therefore see on the pair distribution function a precursor of the interference fringes 
observed when the positions of all the particles are measured! 
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8.1.3 A Monte Carlo simulation 



By sampling the A^— body distribution function P with a Monte Carlo technique, Ja- 
vanainen and Sung Mi Yoo in made a numerical experiment with = 10'^ particles 
and kb = —ka . By distributing the measured positions in a given realization xi,X2, ....xn 



among 30 position bins they obtained histograms like the ones in figure |T6|. It turns out 
that the density in the outcome of each realization of the numerical experiment can be 
fitted by a cosine: 

N 2 

— e^^^^e^^" + e^^^^e^^^ (414) 

where 9a and 9b are phases varying randomly from one realization to the other. In other 
words one has the impression that for each realization the system is in the state 



1 



N 



|0) 



(415) 



with the angle 9 = {9a — 9h)/2 randomly distributed in [— 7r/2,7r/2] . Such a state, 
corresponding to a well defined phase between the two modes a and h , is called a phase 
state 1^ . 



8.1.4 Analytical solution 

We wish to explain the result of the numerical experiment with an analytical argument. 
This has been done with slightly different points of view in |P, We give here what we 
think is the simplest possible presentation. 

Let us allow Poissonian fiuctuations in the number of particles Na and , corre- 
sponding to the distribution probabilities: 



a, b 



(416) 



with mean number of particles Na = Nt = N /2 . These fiuctuations become very small 
as compared to N when the number of particles becomes large: 

AN 1 



for A^, ^ oo. 



The corresponding density operator is a statistical mixture of Fock states: 

oo 

p= Va{Na)Vbm\Na:ka,Nb:kb){Na:ka,Nf,:kb\. 



(417) 



(418) 
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Figure 16: For two different Monte Carlo realizations (a) and (b) of the gedanken experiment, 
histogram of the measured positions of = 1000 particles for an initial Fock state with N/2 
particles in plane wave ka and N/2 particles in plane wave kb = —ka [37|. The positions of 
the particles are expressed in units of 2ii/{ka — kf,) and are considered modulo 2Ti/{ka — kf,) . 



From this form one can imagine that a single realization of the experiment is in a Fock 
state, provided that one keeps in mind that Na and vary in an impredictable way 
from one experimental realization to the other. We known from the work that there 
will be interference fringes in each experimental realization, but this fact is not intuitive. 

The same density operator can also be written in terms of a statistical mixture of 
phase states: 



-N 



d9 



\0)nn{0\. 



(419) 



From this form one can imagine that a single realization of the experiment is in a phase 
state, provided that one keeps in mind that the total number of particles N and the 
relative phase 6 vary in an impredictable way from one realization to the other. This last 
form leads to the following algorithm to generate the positions of the particles according 
to the correct probability distribution: 



1. generate an integer according to the Poisson distribution of parameter N 
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2. generate 6 according to a uniform probability distribution within — 7r/2 and 7r/2 

3. generate the positions Xi, ....xn as if the system was in the state \6)n , in which 
case all the particles are in the same single particle-state and the probability density 
P{xi,....xn) is factorized: 

N 

P{xi,....xn) = l[Pi^j) (420) 

where 

p(x) = ^ le^'^-^e^^ + e^'^^^e-^^p . (421) 

One then obtains interference fringes in each experimental realization, in a very explicit 
way. 

One could also use a third form of the same density operator p , that is a statistical 
mixture of Glauber coherent states: 

p = r — |coh : N:^\^'% coh : iV,^/^e^^^)(coh : iV.^/^e^^% coh : N,'^\^\ 

Jo 2tt Jo 2n 

(422) 

This mathematical form is at the origin of the popular belief that condensates are in 
coherent states. From this form one can only imagine that a single realization of the 
experiment is in a coherent state, keeping in mind that the phases 6a and 6b vary in an 
impredictable way from one realization to the other. In this representation the occurrence 
of interference fringes is straightforward. 

There is an important difference between the coherent states and the Fock or phase 
states: as the number of particles is a conserved quantity in the non-relativistic Hamil- 
tonian used to describe the experiments on atomic gases it seems difficult to produce 
a condensate in a coherent state in some mode ip , that is with p being a pure state 
|coh : a) (coh : a\ where a is a complex number. 

On the contrary one could imagine producing a condensate in a Fock state by mea- 
suring the number of particles in the condensate. One could then obtain a phase state by 
applying a 11/2 Rabi pulse on the Fock state changing the internal atomic state a to a 
superposition (|a) + |6))/v^ where b is another atomic internal state; such a Rabi pulse 
has been demonstrated at JILA and has allowed the measurement of the coherence time 



of the relative phase between the a and b condensates fTQ] • 
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8.1.5 Moral of the story 

• there is in general no unique way of writing the density operator p as a statistical 
mixture. The canonical form corresponding to the diagonalization of p is always 
a possibility but not always the most convenient one. E.g. in our simple model the 
eigenbasis (Fock states) is less convenient than the non-orthogonal family of phase 
states (symmetry breaking states). 

• no measurement or no set of measurements performed on the system can distinguish 
between two different mathematical forms of the same density matrix as a statistical 
mixture. 

• the symmetry breaking point of view consists in writing (usually in an approximate 
way) the A^— body density operator as a statistical mixture of Hartree-Fock states. 
One can then imagine that a given experimental realization of the system is a 
Hartree-Fock state, whose physical properties are immediate to understand as all 
the particles are in the same quantum state. 

• If the system is not in a state that is as simple as a Hartree-Fock state (e.g. in a Fock 
state for our simple model) it is dangerous to make reasonings on the single particle 
density operator (that is on the first order correlation function of the atomic field 
operator) to predict outcomes of single measurements on the system: the relevant 
information may be stored in higher order correlation functions of the field. 



8.2 What is the time evolution of an initial phase state ? 
8.2.1 Physical motivation 

Consider an interference experiment between two condensates A and B either in spa- 
tially separated traps or in different internal states (JlLA-type configuration |7^). Assume 
that the two condensates have been prepared initially in a state with a well defined rel- 
ative phase 6 ; this has actually been achieved at JILA. Let the system evolve freely for 
some time t . How long will the relative phase remain well defined ? This question is 
probably not an easy one to answer. We present here a simple model including only two 
modes of the field. In real life the other modes of the field are not negligible (see for 
example for a discussion of finite temperature effects) and phenomena neglected here 
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such as losses of particles from the trap and fluctuations in the total number of particles 
may be important in a real experiment 0, 

We assume that the state of the system at time t = is a phase state. More 
specifically, expanding the N— th power in Eq. ([415|) with the binomial formula, we take 
as initial state: 

^ / /VI \^/^ 

Mt = 0)) = 2-^/2 Y: WlW^ e^(^-^^)>„ : 0., m : 0.) (423) 

Na=0 \^^a-^^b-/ 

where = N — Na and (f)a,b are the steady state condensate wavefunctions with Na^ 
particles in condensates A, B respectively. The time evolution during t is simple for 
each individual Fock states, as the system is then in a steady state with total energy 

The time evolution of the phase state Eq. ( |423| ) is much more complicated: the state vector 
\^{t)) is a sum of many oscillating functions of time. 

8.2.2 A quadratic approximation for the energy 

The discussion can be greatly simplified if one uses the fact that the binomial factor in 
Eq.( |423| ) for large N is a function of Na and Ni, sharply peaked around Na = Nb = N/2 
with a width y/N : from Stirling's formula n\ ^ (n/e)"\/27m we obtain indeed 

^ ^ ^ ( ^ \ r-^- 1or(NJN)-N, Ior(NJN) ^ ( (N,-Nt)V{2N) 

2^ Njm ^2^\NaNj UN J 

(425) 

We therefore expand the energy E in powers of A^^^ — N/ 2 and Nb — N/ 2 up to second 
order. 

E{Na,Nb) = E{N/2,N/2) + {Na-N/2)dN^E + {Nb-N/2)dN,E 



+i(iV, - N/2Ydl,E + ^{Nb - N/2fdlE 
+ {Na - Nl2){Nb - N/2)dN^dN,E + . . . , (426) 



all the derivatives being taken in {Na,Nb) = {N/2, N/2) . Note that the first derivatives 
of the energy are the chemical potentials fia,b of the two condensates; as the condensates 
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are independent condensates (there is no mechanism locking the relative phase of the 
condensates) one has in general 7^ A^fe • As we restrict to the set of occupation numbers 
such that Na + Nb = N we can rewrite the expansion of the energy using A^^^ — N/2 = 
-{N, - N/2) = {Na - N,)/2 : 

1 fi 

E{Na, N,) ~ E{N/2, N/2) + -(iVa - N,){^^a - fit) + j{Na - N,)\ (427) 
where we have introduced the quantity 

\2 



^ 2h 



^ {dN^-dN,YE . (428) 



8.2.3 State vector at time t 

If one uses the quadratic approximation of the energy the system evolves from the initial 
state Eq. (|423|) to the state 

N / T\T\ \ V2 

\^{t)) = 2-^/2 Y: TFTWl e^(^-^^)(^+''*)e-(^^-^^)^^*/^|iV, : N, : 0,) (429) 

Na=0 V^c^-^^blJ 



The contribution of the term linear in Na — Nf, in Eq. (|427| ) is contained in the quantity 

v = ^{f^b-fia). (430) 

The resulting effect on the time evolution is simply to shift the relative phase between the 
condensate from 6 to 6 + vt : this is a mere phase drift with a velocity v . This phase 
drift takes place only if the 'frequencies' /^a/^ and fib/fi of the atomic fields in A and 
in B are different. 

The effect of the quadratic term in Eq. (|427|) is to spread the relative phase of the 
two condensates. This effect is formalized in ||^, we give here the intuitive result. The 
spreading of the phase can be understood in analogy with the spreading of the wavepacket 
of a fictitious massive particle, with the relative phase 6 being the position x of the 
particle and the occupation number difference Nf, — Na being the wavevector k of the 
particle. The energy term proportional to x plays the role of the kinetic energy of the 
particle responsible for the spreading in position. The effective mass of the fictitious 
particle is M such that 

\(N.~N,f,^§ (431) 
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so that 

M = ^. (432) 

Replacing the discrete sum in Eq. ([429D by an integral we formally obtain the expansion 
of the time dependent state vector of the fictitious particle over the plane waves in free 
space. In this case the variance of the position of the fictitious particle spreads as 

Ax^t) = Ax\0) + (^^J t\ (433) 



Within the approximation (^425|) the wavepacket of the fictitious particle is a Gaussian in 
momentum space, with a standard deviation Ak = (A^/2)^/^ . Initially the position x is 
well defined with a spread ~ 1/ Ak ^ 1 . The relative phase of the condensates will start 
becoming undefined when the position spread Ax of the fictitious particle becomes on 
the order of unity. This happens after a time 

M _ 2^2 

^spread ~ ^ " ^^^2' (^34) 

At times much longer than tgpread it is not correct to replace the discrete sum over 
{Na — Nh)/2 by an integral. The discreteness of A^^ — A^^ leads to reconstructions of a 
phase state (the so-called revivals) at times tg = qn / x , q integer: one can check indeed 
from Eg. ( [429|) that a phase state is reconstructed with a relative phase 6 + vtq + qn /2 for 
A^ even and 9 + vtq for A^ odd. The observability of even the first revival at time ti is 
a non trivial question: the revivals are easily destroyed by decoherence phenomena such 
as the loss of a few particles out of the condensate due to inelastic atomic collisions 0, 
and effects of the non-condensed fraction also need to be investigated. This fragility of 
the revivals is not surprising if one realizes that the state vector |^'(t)) in Eq. ( |429| ) is a 
Schrodinger cat at time ti/2 , that is a coherent superposition of the A^ particles in some 
state 01 and of the A^ particles in some state 02 orthogonal to 0i : the revival at time 
ti is suppressed if the Schrodinger cat at time ti/2 is transformed by decoherence into 
a statistical mixture of the states |A^ : 0i,2) , which is difficult to avoid for large values of 
A^ (see the lecture notes of Michel Brune in this volume). 

8.2.4 An indicator of phase coherence 



To characterize the degree of phase correlation between the two condensates it is natural 
to consider the average of {a)h) where a, h annihilate a particle in condensates A and 
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B respectively. Consider indeed the average over many experimental realizations of some 
one-body observable O sensitive to the relative phase of the two condensates. This 
observable necessarily has a non-vanishing matrix element between the modes (pa and 
(pb so that in second quantized form the part of (O) sensitive to the relative phase involves 
(a'^b) . E.g. in the case of spatially separated condensates one can beat on a 50 — 50 matter 
waves beam splitter atoms leaking out of the condensates and detect the atoms in the 
output channels of the beam splitter 0; the number of counts in the -|- output channel 
averaged of many experimental realizations is proportional to the expectation value of 

O^^^. (435) 

Expanding this product of operators we get 'diagonal' terms such as a^a not sensitive to 
the relative phase, and crossed terms (actually interference terms!) such as sensitive 
to the phase. In the JILA-type configuration, where the condensates A and B are in 
different internal atomic states, an observable O similar to Eq. (|435|) has been achieved 



by mixing the internal states of the two condensates by a 7r/2 electromagnetic pulse and 
by measuring the mean density of atoms in A and B ||70|| . 



From the Schwartz inequality < \ \u\ \ \ \v\\ and setting \u) = a|\E') , \v) = 

we obtain an upper bound for the expectation value of a)b : 

\{^\a^b\^)\ < (^|a^a|^)^/2(^|6^6|^)i/l (436) 

The case of a maximally well defined relative phase corresponds to an equality in this 
inequality, obtained if \u) and \v) are proportional. In the present situation of equal 
mean numbers of particles N/2 in A and in B this corresponds to |\I') being a phase 
state. 

For an initial phase state it is possible to calculate the expectation value of as 
function of time from the expansion (|429|) . One obtains after simple transformations the 
sum 

(atS) it) = ^e--(^+-) E e-^*[^-(^-^« (437) 

with Nf, = N — Na as in Eq . ([429|) . After inspection one realizes that this sum is the 



binomial expansion of a (A^ — 1)— th power so that the final result is 

(at5) (t) = ^e-2*(^+''*) cos^-i xt. (438) 
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From this very simple expression one can calculate the time tc after which the relative 
phase has experienced a significant spread. For short times 1 one can expand the 

cosine function in Eq.( [438| ) to second order in t : 

cos^ Xt = e^'°§™^^* ~ e-^(^*)'/2_ (439) 
One obtains a Gaussian decay of phase coherence with a collapse time 

tc = T-A^ (440) 



\X 



m/2 



equivalent to the rougher estimate Eq. ( [434| ) up to a numerical factor. One can also easily 
see the revivals (reconstruction of |\E') to a phase state) at times tg = qir/x when the 
cosine function is equal to ±1 in Eq.( [438| ). 

Formula ( [44(]| ) can be used to calculate the coherence time of the relative phase of the 
condensates in the present zero-temperature model. As an interesting application of this 
formula we now show that the spreading time of the relative phase can be significantly 
different for mutually interacting and non-mutually interacting condensates. Assume for 
simplicity that the two condensates are stored in cubic boxes of identical size L and with 
periodic boundary conditions. In the MIT-type configuration the two boxes are spatially 
separated and the atoms are in the same internal state; the energy of a configuration with 
Na, Nh atoms in the condensates A, B is then 

E ' 



^Nl + Ni] . (441) 
From Eq. ([428D this form of E leads for an initial phase state to a collapse time of the 



relative phase 

t^ = N'/''— (442) 

where p = N/ {2L^) is the mean spatial density in each of the condensates. In the JILA- 
type configuration the atoms are in the same spatial box but in different internal states; 
the energy of a configuration with Na, Ni, atoms in the condensates A, B is given now 
by Eq.( p94| ) if the two internal states are subject to spatial demixing, or by Eq.( |295| ) if 
there is no demixing instability. The collapse time is then given by 

tc = N^'^- — - (443) 

P[9aa + 9bb - 2[gaa9bbr'^\ 
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for a demixed condensates and by 

tc = N'/^ (444) 

PWaa + Qbh - '2^9 ab\ 

for fully overlapping condensates. When the coupling constants among the various internal 
states are close to each other the denominator in Eqs. ( [443| , [444| ) can become small, which 
results in a relative phase coherence time tc much larger than in the MIT-configuration 
Eg ■ ( |442| ) ■ This fortunate feature of close coupling constants is present for rubidium in the 
JILA experiment | |7D| ! 

In real life the condensates are usually stored in harmonic traps; the simple formulas 
obtained for a cubic box have to be revisited. This has been done analytically for spatially 
separated condensates and numerically for mutually interacting condensates 

m. 



9 Symmetry breaking description of condensates 

We have already seen in chapter | that it is very convenient, physically, to introduce 
phase states to understand the phenomenon of interference between two Bose-Einstein 
condensates: rather than assuming that two Bose-Einstein condensates that "have never 
seen each other" are in Fock states, one assumes that they are in a phase state with a 
relative phase varying in an unpredictable way for any new experimental realization. One 
can even suppose that the condensates are in coherent states of the atomic field; this 
description is said to 'break the symmetry', here the U{1) symmetry associated to the 
invariance of the Hamiltonian by a change of the phase of the atomic field operator. 

In this chapter we consider other examples of symmetry breaking descriptions: 5*0(3) 
symmetry breaking (case of spinor condensates) and spatial translational symmetry break- 
ing (case of one dimensional condensates with attractive interactions). In both cases the 
procedure is the same: the ground state of the system is symmetric, its mean-field ap- 
proximation by Hartree-Fock states breaks the symmetry. In both cases we will consider 
Gedanken experiments whose single outcomes can be predicted from the exact ground 
state and from the Hartree-Fock state. This will illustrate the ability of the mean-field 
approximation to allow physical predictions in an easy and transparent way, correct in 
the limit of a large number of particles. 
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9.1 The ground state of spinor condensates 

The alkali atoms used in the Bose-Einstein condensation experiments have an hyperfine 
structure in the ground state, each hyperfine level having several Zeeman sublevels. We 
have up to now ignored this structure in the lecture, as we were implicitly assuming that 
the atoms were polarized in a well defined Zeeman sublevel. 

Consider for example Na atoms used at MIT in the group of Wolfgang Ketterle. 
The ground state has an hyperfine splitting between the lower multiplicity of angular 
momentum F = 1 and the higher multiplicity of angular momentum F = 2 . All the 
three Zeeman sublevels mp = 0, ±1 of the lower multiplicity F = 1 cannot be trapped 
in a magnetic trap (if mp = —1 is trapped than mp = +1 which experiences an 
opposite Zeeman shift is antitrapped). But they can all be trapped in an optical dipole 
trap, produced with a far off-resonance laser beam, as the Zeeman sublevels experience 
then all the same lightshift. This optical trapping was performed at MIT [^, opening 
the way to a series of interesting experiments with condensates of particles of spin one 



We concentrate here on a specific aspect, the ground state of the spinor condensate, 
assuming for simplicity that the atoms are stored in a cubic box with periodic boundary 
conditions. 

9.1.1 A model interaction potential 

We have to generalize the model scalar pseudo-potential of Eq. (|73D to the case of particles 
having a spin different from zero. As we want to keep the simplicity of a contact interaction 
potential we choose the simple form 



that is the product of an operator acting only on the spin of the particles 1 and 2, and 
of the usual regularized contact interaction acting only on the relative motion of the two 
particles. The interaction potential V{1, 2) has to be invariant by a simultaneous rotation 
of the spin variables and of the position variables of the two particles. As the contact 
interaction is already rotationally invariant, the spin part of the interaction Vspin(l,2) 
has to be invariant by any simultaneous rotation of the two spins. 



77 . 




(445) 
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This condition of rotational invariance of Vspin(l, 2) is easy to express in the coupled 
basis obtained by the addition of the two spins of particle 1 and particle 2: within each 
subspace of well defined total angular momentum Vspin(l,2) has to be a scalar. Let us 
restrict to the case studied at MIT, with spin one particles. By addition of F = 1 and 
F = 1 we obtain a total angular momentum Fi^t = 2 , 1 or 0, so that one can write 

Vspin(l, 2) = ^2Pnot=2(l, 2) + ^i^not=i(l' 2) + goPF,,,=o{l, 2) (446) 

where the g 's are coupling constants and the P(l,2) 's are projectors on the subspace 
of particles 1 and 2 with a well defined total angular momentum Ftot • At this stage we 
can play a little trick, using the fact that the states of Ftot = 1 are antisymmetric by the 
exchange of particles 1 and 2 (whereas the other subspaces are symmetric). The regu- 
larized contact interaction scatters only in the s -wave, where the external wavefunction 
of atoms 1 and 2 is even by the exchange of the positions fi and r2 ; as our atoms are 
bosons, the spin part has also to be symmetric by exchange of the spins of atoms 1 and 2 
so that the 'fermionic' part of Vspin(l, 2) , that is in the subspace Ftot = 1 , has no effect. 
We can therefore change gi at will without affecting the interactions between bosons. 
The most convenient choice is to set gi = g2 so that we obtain 

Vspi„(l, 2) = ^2ld(l, 2) + (^0 - ^2)Pnot=o(l> 2) (447) 

where Id is the identity. The subspace Ft^t = is actually of dimension one, and it is 
spanned by the vanishing total angular momentum state |-?/'o(li2)) . Using the standard 
basis |m = —1, 0, +1) of single particle angular momentum with z as quantization axis, 
one can write 

|^o(l, 2)) = -i= [| + 1, -1) + I - 1, +1) - |0, 0)] . (448) 

A more symmetric writing is obtained in the single particle angular momentum basis 
\x, y, z) used in chemistry, defined by 

1 + 1) = -l=^\x)+t\y)) (449) 

1-1) = +l={\x)-t\y)) (450) 
|0) = 1^). (451) 
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The vector |a) in this basis (a = x, y, z) is an eigenvector of angular momentum along 
axis a with the eigenvalue zero. One then obtains 

iV'oll, 2)) = [Ix, x) + \y, y) + \z, z)] . (452) 

To summarize the part of the Hamiltonian describing the interactions between the 
particles can be written, if one forgets for simplicity the regularizing operator in the 
pseudo-potential: 



^int = y / E V'aV'JV'/3V'a 

+ ^^/rf'r i^i^ii^P^P- (453) 

where V'a(^) is the atomic field operator for the spin state \(y) ■ This model Hamiltonian 
has also been proposed by |]7S|, |7S|, PU[ . 

9.1.2 Ground state in the Hartree-Fock approximation 

As we are mainly interested in the spin contribution to the energy we assume for simplicity 
that the condensate is in a cubic box of size L with periodic boundary conditions. We 
assume that the interactions between the atoms are repulsive {g-2,gQ>Q) and we suppose 
that there is no magnetic field applied to the sample. 

We now minimize the energy of the condensate within the Hartree-Fock trial stat- 
evectors \Nq : 0) with the constraint that the number of particles A^o is fixed ( |0) is 
normalized to unity) but without any constraint on the total angular momentum of the 
spins. The external part of the condensate wavefunction is simply the plane wave with 
momentum = whereas the spinor part of the wavefunction remains to be determined: 

^''\^) = Tm E with 5]|c„|2 = l. (454) 

a=x,y,z a 

Using the model interaction Hamiltonian Eq. ( |453| ) we find for the mean energy per particle 
in the condensate 

92 + ^rp^igo - 92)\A\ (455) 



No 2L3 6L3 
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where we have introduced the complex quantity 

A= Y: cl = c' (456) 

a=x,y,z 

where c is the vector of components {c^, Cy,Cz) ■ We have to minimize the mean energy 
over the state of the spinor. 

• Case g2 > go 



This is the case of sodium |j7^. As the coefficient qq — g2 is negative in Eg .( [4551 ) we have 



to maximize the modulus of the complex quantity A . As the modulus of a sum is less 
than the sum of the moduli we immediately get the upper bound 

1^1 < E Ica|' = l (457) 

a=x,y,z 

leading to the minimal energy per particle 

J^^=^^92 + ^^{%-9.). (458) 

The upper bound for \A\ is reached only if all complex numbers c\ have the same phase 
modulo 27r . This means that one can write 

Ca = e^^n„ (459) 

where 9 is a constant phase and n = {nx,ny,nz) is any unit vector with real components. 
Physically this corresponds to a spinor condensate wavefunction being the zero angular 
momentum state for a quantization axis pointing in the direction of n . The direction 
of n is well defined in the Hartree-Fock ansatz, but it is arbitrary as no spin direction 
is privileged by the Hamiltonian. We are facing symmetry breaking, here a rotational 
S0{3) symmetry breaking, as we shall see. 

• Case g2 < go 

In this case we have to minimize \A\ to get the minimum of energy. The minimal value 
of \A\ is simply zero, corresponding to spin configurations such that 

c'^ E 4 = (460) 

a=x,y,z 
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with an energy per condensate particle 

E No -I 



92. (461) 



No 2L3 

To get more physical understanding we split the vector c as 

c = R + iI (462) 

where the vectors R and / have purely real components. Expressing the fact that 
the real part and imaginary part of vanish, and using the normalization condition 
c ■ c* = 1 in Eg. ( |454| ) we finally obtain 



R-I = (463) 
R^ = I'^ = \- (464) 

This means that the complex vector c is circularly polarized with respect to the axis Z 
orthogonal to I and R . Physically this corresponds to a spinor condensate wavefunction 
having an angular momentum ±^ along the axis Z . The direction of axis Z is well 
defined in the Hartree-Fock ansatz but it is arbitrary. 

9.1.3 Exact ground state of the spinor part of the problem 

Imagine that we perform some intermediate approximation, assuming that the particles 
are all in the ground state = of the box but not assuming that they are all in the 
same spin state. We then have to diagonalize the model Hamiltonian 

+ 7rT^Ago-92)A^A (465) 

a,l3=x,y,z 

where da annihilates a particle in state |A; = 0)|q;) { a = x,y, z ) and where we have 
introduced 

A = dl + dl + dl (466) 

Up to a numerical factor A annihilates a pair of particles in the two-particle spin state 
|'?/'o(l, 2)) of vanishing total angular momentum, as shown by Eq. ([452|) . 

The Hamiltonian Eq.( [465| ) can be diagonalized exactly |8l|. This is not surprising as 



(i) it is rotationally invariant and (ii) the bosonic A^o~ particle states with a well defined 
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total angular momentum S^o can be calculated: one finds that S^^ = N^, Nq — 2, . . . , 
leading to degenerate multiplicities of ifspin of degeneracy 2Sno + ^ • 

In practice one may use the following tricks: The double sum proportional to g2 in 
Eg .( [4651 ) can be expressed in terms of the operator number Nq of condensate particles 
only, 

^o = E4«a. (467) 

a 

So diagonalizing -ffspin amounts to diagonalizing AM ! 

Second the total momentum operator S of the Nq spins, defined as the sum of all 
the spin operators of the individual atoms in units of h , can be checked to satisfy the 
identity 

S-S + A^A = No{No + l) (468) 



so that the Hamiltonian for A^^o particles becomes a function of 5* |81 



Or, ^ ^ 1 



Nq{No + ^)-S-S 



(469) 



We recall that S ■ S = S^fX^No + 1) within the subspace of total spin S^g . 

When g2 < go the ground state of -ffspin corresponds to the multiplicity Snq = , 
containing e.g. the state with all the spins in the state |+) . In this case the A^o~ particle 
states obtained with the Hartree-Fock approximation are exact eigenstates of -f^spin • 

When g2 > go the ground state of ifspin corresponds to the multiplicity of minimal 
total angular momentum, Snq = 1 for No odd or Snq = for A^o even. In this case 
the Hartree-Fock state is a symmetry breaking approximation of the exact ground state 
of -ffspin • The error on the energy per particle tends to zero in the thermodynamical 
limit; for A^^o even one finds indeed 

SE 1 

^ = -3jjtoo-9.). (470) 

But what happens if one restores the broken symmetry by summing up the Hartree- 
Fock ansatz over the direction n defined in Eg . (14591) ? Assume that A'q is even; one 
has then to reconstruct from the Hartree-Fock ansatz a rotationally invariant state. This 
amounts to considering the following normalized state for the A"o spins: 

72 



\M,) = ^No + l I ^\No:n) (471) 
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where (Pn indicates the integration over the unit sphere (that is over all solid angles) 
and \No : fi) is the state with Nq particles in the single particle state 

\n) = nr,\x) +ny\y) +n^\z). (472) 

The state vector |\E') , being non zero and having a vanishing total angular momentum, 
is equal to the exact ground state of i^spin ! 

The expression (|471|) can be used as a starting point to obtain various forms of |\E') . 
If one expresses the Hartree-Fock state as the Nq -th power of the creation operator 
J2a ^I'^a acting on the vacuum |vac) , and if one expands this power with the usual 
binomial formula, the integral over n can be calculated explicitly term by term and one 
obtains: 

1^) =Ar(it)'^°/'|vac) (473) 

where A/" is a normalization factor and the operator A is defined in Eq.( [466| ). Formula 
(|473|) indicates that l^') is simply a 'condensate' of pairs in the pair state \ipo{'^, 2)) • It 
can be used to expand |\E') over Fock states with a well defined number of particles in 
the modes m = 0,m = ±1 , reproducing Eq.(13) of |PT| . 

To be complete we mention another way of constructing the exact eigenvectors and 
energy spectrum of i/spin • The idea is to diagonalize A^A using the fact that A obeys 
a commutation relation that is reminiscent of that of an annihilation operator: 

[A, A^] = ANo + 6. (474) 

In this way A"^ acts as a raising operator: acting on an eigenstate of A'^A with eigenvalue 
A and Nq particles, it gives an eigenstate of A^A with eigenvalue X + ANq + G and with 
A^o + 2 particles. One can also check from the identity ( [46 8|) that the action of A^ does 
not change the total spin: 

[A\S-S] = 0. (475) 

By repeated actions of A^ starting from the vacuum one arrives at Eq.(^73|), creating 
the eigenstates with A^^o even and vanishing total spin S = . By repeated actions of 
A^ starting from the eigenstates with A^^q = 2 and total spin S = 2 (e.g. the state 
I + +) ) one obtains all the states with A^q even and total spin S = 2 . More generally 
the eigenstate of i^spin with total spin S , a spin component m = S along z and A'o 
particles is: 

\\NQ,S,m = S) oc n^y^''-^^^^ \s : +1) (476) 
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where IS* : +1) represents S particles in the state | + 1) . From Eg. ( |476| ) one can 
generate the states with spin components m = S — 1, . . . , —S by repeated actions of the 
spin-lowering operator S- = — iSy in the usual way. We note that formula ([476|) was 
derived independently in [Q. 

9.1.4 Advantage of a symmetry breaking description 

Imagine that we have prepared a condensate of sodium atoms { g2 > go ) m the collective 
ground spin state, and that we let the atoms leak one by one out of the trap, in a way 
that does not perturb their spin. We then measure the spin component along z of the 
outgoing atoms. Suppose that we have performed this measurement on k atoms, with 
k <ti No . We then raise the simple question: what is the probability pk that all the k 
detections give a vanishing angular momentum along z ? 

Let us start with a naive reasoning based on the one-body density matrix of the 
condensate (even if the reader has been warned already in § 8.1.2| on the dangers of such 



an approach!). The mean occupation numbers of the single particle spin states |m = — 1) , 
|m = 0) and |m = +1) in the initial condensate are obviously all equal to iVo/3 , as the 
condensate is initially in a rotationally symmetric state. The probability of detecting the 
first leaking atom in |m = 0) is therefore 1/3 . Naively we assume that since k Nq 
the detections have a very weak effect on the state of the condensate and the probability 
of detecting the n -th atom ( n < ) in the m = channel is nearly independent of the 
77, — 1 previous detection results. The probability for k detections in the m = channel 
should then be 

Pr^^ = ^- (477) 

Actually this naive reasoning is wrong (and by far) as soon as k > 2 . The first 
detection of an atom in the m = channel projects the spin state of the remaining 
atoms in 

l^i) = Alaol^) (478) 

where o-o annihilates an atom in spin state m = , is the collective spin ground 
state (|7T) and Afi is a normalization factor. The probability of detecting the second 



atom in m = (knowing that the first atom was detected in m = ) is then given by 
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The denominator is simply equal to A^^q — 1 as |\E'i) is a state with Nq — 1 particles. Using 
the integral form ( |471|) and the simple effect of an annihilation operator on a Hartree-Fock 
state, e.g. 

al\No : n) = [No{No - l)]'/'n^|iVo - 2 : n) (480) 
we are able to express the probability in terms of integrals over solid angles: 

= ^-r-^r (481) 

Pi d^n d^n' n,n',{n-n'f°'' 

We suggest the following procedure to calculate these integrals. One first integrates over 
n' for a fixed n , using spherical coordinates relative to the 'vertical' axis directed along 
n : the polar angle 6' is then the angle between n' and n so that one has simply 
n- n' = cos 6' . The integral over 6' and over the azimuthal angle 0' can be performed, 
giving a result involving only . The remaining integral over n is performed with the 
spherical coordinates of vertical axis z . This leads to 



The ratio P2/P1 is therefore different from the naive (and wrong!) prediction (477) 



For A^o = 2 one finds = 1 so that the second atom is surely in m = if the first 

atom was detected in m = . As the two atoms were initially in the state with total 
angular momentum zero, this result could be expected from the expression ( [448|) of the 



two-particle spin state. In the limit of large Nq we find that once the first atom has been 
detected in the m = channel, the probability for detecting the second atom in the same 
channel m = is 3/5 . This somehow counter-intuitive result shows that the successive 



detection probabilities are strongly correlated in the case of the spin state (|471|) 
The exact calculation of the ratio 



Pk+i 



d^n I d^n' n'l+^n'^+Hn ■ nY"-^''^^^ 

(483) 



Pk d^n d''n'nln'^{n-n'f°-'' 



is getting more difficult when k increases. The large A^o limit for a fixed k is easier to 
obtain: in the integral over n' the function (n- n')^"~('^+^) is extremely peaked around 
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n' = n so that we can replace n'^^^ by n^^^ . This leads to 

hm ^ = ^. (484) 

We now give the reasoning in the symmetry breaking point of view, which assumes 
that a single experimental realization of the condensate corresponds to a Hartree-Fock 
state lA^^o : n) with the direction n being an impredictable random variable with uniform 
distribution over the sphere. If the system is initially in the spin state |A^o '■ ^) there 
is no correlation between the spins, and the probability of having k detections in the 
channel m = is simply {n1)^ . One has to average over the unknown direction n to 
obtain ^ ^ 



pf = / :_nf = -. (485) 

One recovers in an easy calculation the large Nq limit of the exact result. Eg. ( |484]) ! We 
note that the result ( [48 5| ) is much larger than the naive (and wrong) result (|477|) as soon 
as /c ^ 1 . 



9.2 Solitonic condensates 

We consider in this section a Bose-Einstein condensate with effective attractive inter- 
actions subject to a strong confinement in the x — y plane so that it constitutes an 
approximate one- dimensional interacting Bose gas along z . Such a situation is interest- 
ing physically as it gives rise in free space to the formation of 'bright' solitons well known 
in optics but not yet observed with atoms. Also the model of a one-dimensional Bose gas 
with a 6 interaction potential has known exact solutions in free space, that can be used 
to test the translational symmetry breaking Hartree-Fock approximation. 



9.2.1 How to make a solitonic condensate ? 

Consider a steady state condensate with effective attractive interactions in a three di- 
mensional harmonic trap. The confinement in the x — y plane is such that the trans- 
verse quanta of oscillation hu^^y are much larger than the typical mean field energy per 
particle A'"o|(7||0p, where is the condensate wavefunction with Nq particles. This 
confinement prevents the occurrence of a spatial collapse of the condensate (see § ^.2.1|) . 



Broken symmetry 



131 



The confinement is however not strong enough to violate the validity condition of the 
Born approximation for the pseudo-potential, k\a\ <C 1 with k ~ [muJx^y/fiY^'^ ■ 

In this case we face a quasi one- dimensional situation, where the condensate wave- 
function is approximately factorized as 



(486) 



where Xx and Xy a-^'s the normalized ground states of the harmonic oscillator along x 
and along y respectively. By inserting the factorized form (|486|) in the Gross-Pitaevskii 



energy functional Eq. (|139|) and by integrating over the directions x and y we obtain an 
energy functional for ijj : 



No I dz 



2m 



dip 



dz 



+ ^mLuy\^iz)\' + ^NogMm^)\' 



(487) 



where we have dropped the zero-point energy of the transverse motion and we have called 
gid the quantity 



gid = g dx dy \Xx{x)\ IXyivW 



9 



2'Kh 



The corresponding time independent Gross-Pitaevskii equation for is 



iii){z) 



f^dSp_ 

2m dz"^ 



+ 



-muy + NogM\^{z)\^ 



(488) 



(489) 



The energy functional Eq.( [487| ) corresponds to a one-dimensional interacting Bose gas 
with an effective coupling constant between the atoms equal to gid , that is one can 
imagine that the particles have a binary contact interaction 



V{zi,z2) = gid^zi - Z2) 



(490) 



Note that such a Dirac interaction potential leads to a perfectly well defined scattering 
problem in one dimension, contrarily to the three dimensional case. 

Imagine now that we slowly decrease the trap frequency along z while keeping intact 
the transverse trap frequencies, until Uz vanishes. What will happen then? If g was pos- 
itive the cloud would simply expand without limit along z . With attractive interaction 
the situation is dramatically different: due to the slow evolution of uj^ the condensate 
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wavefunction will follow adiabatically the minimal energy solution of the Gross-Pitaevskii 
equation. For Uz = this minimal energy solution is the so-called bright soliton, well 
known in non-linear optics. We recall the analytic form of the solitonic wavefunction: 

'"<^) = MVS js;^ '''''' 



where / is the spatial radius of the soliton: 



/ = • (492) 

Nomgid 

Note that this size / results of a compromise between minimization of kinetic energy by 
an increase of the size and minimization of interaction energy by a decrease of the size, so 
that the typical kinetic energy per particle h'^/ {mt^) is roughly opposite to the interaction 
energy per particle Nogid/l . We also give the corresponding chemical potential: 



We briefly address the validity of the Gross-Pitaevskii solution ( [49 1| ). As we have 
pointed out in the three dimensional case (see for example § |3.2.1| ) we wish that the Born 
approximation for the interaction potential be valid. In one dimension the 6 interaction 
potential can be treated in the Born approximation only if the relative wavevector of the 
colliding particles is high enough (in contrast to the three-dimensional case): 



mgid 



> 1. (494) 



This condition can be obtained of course from a direct calculation, but also from a dimen- 
sionality argument ( mgid/h^ is the inverse of a length) and from the fact that the Born 
approximation should apply in the limit gid for a fixed k . If we use the estimate 
^ 1// we obtain the condition 

~ iVo > 1, 495 

mgiJ 

implicitly valid here as we started from a condensate! 

Another phenomenon neglected in the prediction (|491|) is the spreading of the center 
of mass coordinate during the switch-off of the trapping potential along z . Whereas 
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Eq.( [491| ) assumes that the abscissa of the center of the sohton zq is exactly the 
spreading of the center of mass leads in real life to a finite width probability distribution 
for zo . This spreading can be calculated simply for an almost pure condensate iVo ~ iV , 
using the fact that the center of mass coordinate operator Z and the total momentum 
operator P of the gas along z axis are decoupled from the relative coordinates of the 
particles in a harmonic potential, in presence of interactions depending only on the relative 
coordinates. To prove this assertion one expresses the operators Z and P in terms of 
the position and momentum operators of each particle i of the gas: 

1 ^ 

Z = j^^z, (496) 

N 

P = T.P^ (497) 
1=1 

and one derives the following equations of motion in Heisenberg point of view: 

dZ P , 

(498) 



dt Nm 
dp 

^ = -Nmujl{t)Z. (499) 

The spreading acquired by Z is not negligible when it becomes comparable to the size / 
of the soliton. 

The spreading of Z is interesting to calculate in the absence of harmonic confinement 
along z , ujz = , with the simple assumption that all the particles of the gas are at time 
t = in the soliton state IV') of Eq.( |491| ). As P is a constant of motion for Uz = 
one has simply 

Pt 

Zit) = ZiO) + ^ (500) 
so that the variance of the center of mass coordinate at time t is 

Var(Z)(t) = Var(Z)(0) + ^(^(0)P + P^(0)) + ^Var(P). (501) 



One then replaces Z{0) and P by the sums ( [496| , [49 7|) . As the single particle wave- 



function ip has vanishing mean position and mean momentum all the 'crossed terms' 
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expectation values involving two different particles vanish. As iplz) is a real wavefunc- 
tion one finds also {tp\zp + pz\tp) = so that the contribution linear in time vanishes. 
One is left with 

" ^\mp'\^). (502) 



Var(Z)(t) = -(^|.^|^) + ^^^ 

The variance of Z , initially times smaller than the single particle variance {iplz"^ 
becomes equal to the single particle variance after a time 



tc 



1/2 



2n 



(503) 



where we used the explicit expressions 



2;2 



12 



3/2' 



(504) 
(505) 



The spreading phenomenon of the position of the soliton is formally equivalent to the 
spreading of the relative phase of two condensates initially prepared in a phase state (see 
The critical time tc in ( |503|) scales as N'^/'^fi/\iJi\ as in Eq. ([442|) . 



9.2.2 Ground state of the one-dimensional attractive Bose gas 

We consider here the model of the one-dimensional gas of bosonic particles interacting 
with the contact potential Eq . (|490|) and in the absence of any confining potential. 

It turns out that in this model with gid > one can calculate exactly the eigenenergies 
and eigenstates of the Hamiltonian for A^ particles using the Bethe ansatz |8^. We 
consider here the less studied attractive case gid < , where several exact results are also 
available. In particular the exact expression for the ground state energy is known MM: 



EoiN) 



1 mgld 

'24 



I] 



(506) 



and the corresponding A^— particle wavefunction of the ground state is [B5 



^(^1, 



Af exp 



mgid 
2h^ 



E 



(507) 



l<i<j<N 
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To determine the normalization factor Af we enclose the gas in a fictitious box of size L 
tending to +cxd : ^ 

,2 (iV-1)! (m\gu\Y'' 



(508) 



NL \ j 

To what extent can we recover these results using a Hartree-Fock ansatz : ^z^) for 
the ground state wavefunction? As discussed around Eq. (|151|) we get a mean energy for 
the Hartree-Fock state very similar to Eq.( |487| ): 



N / dz 



2m 



dip 



dz 



+ -{N - l)g,,mz)\' 



(509) 



We minimize this functional using the results of §9.2.1| , replacing Nq by — 1 , and we 



obtain 



1 mgja 



(510) 



24 

The deviation of the Hartree-Fock result from the exact result is a fraction of the 

energy and is small indeed in the large N limit, as expected from the validity condition 
(|495|) ! 

There is a notable difference of translational properties however. Whereas the exact 
ground state (|507|) is invariant by a global translation of the positions of the particles, 
as it should be, the Hartree-Fock ansatz leads to condensate wavefunctions ip localized 
within the length / around some arbitrary point Zq (around = in Eq.( [491| )): 

1 1 



(2/)V2 cosh[{z- Zo)/l] 



(511) 



2r 



(512) 



with a spatial radius 

^ = -7 

(A^ - l)mgM 

The Hartree-Fock ansatz \N : ip) therefore breaks the translational symmetry of the 
system. 

^ The center of mass of the gas corresponds to a fictitious particle of wavevector K , where hK 
is the total momentum of the gas, and of position Z , where Z is the centroid of the gas. In the 
ground state \'^) the center of mass is completely delocalized with K ^ . The factor 1/L in jA/'P 
originates from the normalization of the fictitious particle plane wave in the fictitious box of size L , 



{Z\K) 



iKZ 



I \/L . The more correct mathematical way (not used here) is to normalize in free space 



(no box) using the closure relation J dK\K){K\ = Id , which amounts to replace L by 27r 
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Breaking a symmetry of the system costs energy, and this can be checked for the 
present translational symmetry breaking. As the center of mass coordinates Z, P of the 

particles are decoupled from the relative coordinates of the particles we can write the 
total energy of the gas as the sum of the kinetic energy of the center of mass and an 
'internal' energy including the kinetic energy of the relative motion of the particles and 
the interaction energy. Whereas the exact ground state wavefunction has a vanishing 
center of mass kinetic energy, the symmetry breaking ansatz \N : ip) contains a center 
of mass kinetic energy: 

E,,,,^, = {N : ^\^\N : i;) (513) 

where mN is the total mass of the gas and P is the total momentum operator. Using 
the definition ( [4971 ), expanding the square of P , and using the fact that the soliton 
wavefunction ip has a vanishing mean momentum we obtain 

2 

E,.o.m. = {^\^\^) (514) 

= ^^(^-1)'- (515) 

We see that -Ec.o.m. accounts for half the energy difference between the exact ground state 
energy ( |506| ) and the Hartree-Fock energy (|510|) . 



9.2.3 Physical advantage of the symmetry breaking description 

We now raise the question: is there a Bose-Einstein condensate in the one-dimensional 
free Bose gas with attractive interaction? To make things simple we assume that the 
gas is at zero temperature so that the A^— particle wavefunction is known exactly, see 
Eq.(pO^). 



We start with a reasoning in terms of the one-body density operator (even if we know 
from the previous physical examples that this may be dangerous). Paraphrasing the 
usual three dimensional definition of a Bose-Einstein condensate in free space we put the 
one-dimensional gas in a fictitious box of size L and we calculate the mean number of 
particles Uq in the plane wave with vanishing momentum p = in the limit L +oo . 

The calculation with the exact ground state wavefunction has been done [|^. One 
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finds that no is going to zero as 1/L 



2h^ 



lim noL = C(N)— r. (516) 

.^+00 rn\gid\ 



The factor C{N) is given by 



N N 

I '/ I II I /V 7 I 1 I 

(517) 



{j -ly. {N ^ [ 1 



k{N+l-k)--{N + l) 



and converges to 7r^/2 in the large hmit, so that Uq no longer depends on in 
this limit. There is therefore no macroscopic population in the p = momentum state. 
One may then be tempted to conclude that there is no Bose-Einstein condensate, even at 
zero temperature, in the one-dimensional Bose gas with attractive contact interactions. 
However we have learned that a reasoning based on the one-body density matrix may 
miss crucial correlations between the particles, and that the symmetry breaking point of 
view may be illuminating in this respect. 

The translational symmetry breaking point of view approximates the state of the gas 
by the -body density operator: 

rL/2 Ap'n 

p^^= \im / ^\N:i;J{N:^,,\. (518) 



J-L/2 L 

In the large limit we expect this prescription to be valid for few-body observables. Of 
course for a A^— body observable such as the kinetic energy of the center of mass of the 
gas, the results will be different, Eq.( |515|) for the symmetry breaking point of view vs. a 
vanishing value for the exact result. 

Let us test this expectation by calculating in the Hartree-Fock approximation the mean 
number of particles in the plane wave {z\k) = exp(ikz) / L^^"^ . Using the following action 
of the annihilation operator of a particle with wavevector k on the Hartree-Fock 
state: 

ak\N : = N'/^k\^J\N - 1 : (519) 

we obtain 

nl' = N\{km'. (520) 

The momentum distribution of the particles in the gas in this approximation is simply 
proportional to the momentum distribution of a single particle in the solitonic wavefunc- 
tion ip ! It turns out that the Fourier transform of the 1 / cosh function can be calculated 
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exactly, and it is also a 1/ cosh function. We finally obtain: 

hf 1 ^^^^ 1 

where / is the soliton size given in Eg .( ^121) . For k = one recovers the large N limit 
of the exact result 



In more physical terms, one can imagine from Eq.(518) that a given experimental 
realization of the Bose gas corresponds to a condensate of particles in the solitonic 
wavefunction ( |511| ), with a central position zq being a random variable varying in an 
unpredictable way for any new realization of the experiment. There is therefore a Bose- 
Einstein condensate in the one- dimensional attractive Bose gas! 

An illustrative gedanken experiment would be to measure the positions along z of all 
the particles of the gas. In the symmetry breaking point of view the positions zi, . . . ,zn 
obtained in a single measurement are randomly distributed according to the density 
= \iIj{z — zq)\'^ where zq varies from shot to shot as the relative phase of the two 
condensates did in the MIT interference experiment. As we know the exact ground state 
(|507|) we also know the exact A^— body distribution function, |^E'(2;i, . . . , Z]sf)\'^ . This is 
however not so easy to use! 

So we suggest instead to consider the mean spatial density of the particles knowing 
that the center of mass of the cloud has a position Z . In the exact formalism this gives 
1: 



p{z\Z) = /cizi... /(i2^|^(zi,...,zjv)n 





where / is the -dependent length of the soliton ( pl2| ), the integrals are taken in 
the range [—L/2, L/2] and L +oo ; the factor L , compensating the one in the 
normalization factor of , ensures that the integral of p{z\Z) over z is equal to . 

In the symmetry breaking point of view the definition of p{z\Z) is similar to Eq. (|522|) ; 
the factor L cancels with the 1/L factor of Eq. (|518|) . This leads to 



1 ^ 



p'\z\Z) = JdzoJdz,..JdzN(j[mzk-zo)\''^(^5iz-z,)^5^^^ 



n=l 
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N J dz^...j dzN[j{mzk)\^^5l^Z-z + z^-^f^^z}i (523) 



where we have made the change of variables zi^ + zq (which allows to integrate 

over Zq ) and we have replaced the sum over the indiscernible particles j by times 
the contribution of particle j = 1 . The multiple integral over the positions zi,...,zi^ 
can be turned into a single integral over a wavevector q by using the identity ^{X) = 
J dq/{2n) exp(igX) , allowing a numerical calculation of p^'^(2;|Z) . 

Does the approximate result ( |523|) get close to the exact result for large ? We 
compare numerically in figure ^ the exact density p{z\Z) to the symmetry breaking 
mean-field prediction p^^{z\Z) : modestly large values of give already good agree- 
ment between the two densities. This validates the symmetry breaking approach for the 
considered gedanken experiment. 

What happens in the large N limit? In Eg. ( p23D each variable explores an interval 
of size ~ / so that the quantity {zi + . . . + Z]\f)/N has a standard deviation ~ l/y/N 
much smaller than I and can be neglected as compared to zi inside the 6 distribution. 
This leads to 

p'^^izlZ) - N\^,,=z{z)\^ for yiV>l (524) 

where the solitonic wavefunction ipzo=z is given in Eg. ( [51 1] ). Numerical calculation of 
p^^{z\Z) shows that Eg. ( |524| ) is a good approximation over the range \z — Z\ ~ / for 
N = 10 already! 
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Figure 17: For the ground state of the one-dimensional attractive Bose gas, position dependence 
of the mean density of particles knowing that the center-of-mass of the gas is in Z = . Solid 
line: exact result p(z\Z = 0) . Dashed line: mean-field approximation p^^{z\Z = 0) . The 
position z is expressed in units of the 'soliton' radius / given in Eq.( |512| ), and the linear 
density in units of N/l . The number of particles is (a) = 10 and (b) = 45 . 
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